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Abstract 

Collisional  cross  sections  were  calculated  as  a  function  of  energy  for  two  coupled  one 
dimensional,  spherically  symmetric  potentials.  The  Split  Operator  Method[4]  was 
used  to  propagate  an  initial  Mpller  state,  chosen  to  be  a  Gaussian  in  the  asymptotic 
limit,  through  a  potential.  The  correlation  between  the  wave  packet  and  Mpller  hnal 
state  was  calculated  at  each  time  step.  Using  the  Ghannel  Packet  Method[18],  the 
correlation  function  was  used  to  obtain  scattering  matrix  elements.  From  scattering 
matrix  elements  for  several  different  effective  potential  values  and  using  the  Method  of 
Partial  Waves  [7],  the  collisional  cross  section  is  calculated  for  the  transition  from  the 
^P3/2  to  ^Pi/2  level.  This  method  was  applied  to  LiHe  at  low  energy,  with  results  close 
to  experimentally  measured  values  for  Alkali-Noble  Gas  interactions.  Gross  sections 
were  also  calculated  for  NaHe,  KHe,  RbAr,  and  LiHe-3.  An  explanation  for  the  low 
RbAr  cross  sections  from  the  time  dependent  calculations  is  provided. 
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CALCULATION  OF  COLLISIONAL  CROSS  SECTIONS  FOR  THE 
2p3/2  ^  '^Pxi2  TRANSITION  IN  ALKALI-NOBLE  GAS  SYSTEMS 

I.  Overview 


1.1  Introduction 

Collisional  cross  sections  play  a  fundamental  role  in  atomic  and  molecular  transi¬ 
tions  in  many  systems.  They  provide  a  non-radiative  source  of  energy  transfer  that 
is  particularly  important  in  laser  systems  [20] .  DPALs  are  one  area  of  interest  in  the 
Air  Force  where  collisions  are  essential  to  operation  of  the  laser  [14] .  The  energy  level 
diagrams  for  all  alkali  atoms  are  similar  in  structure,  although  they  differ  in  exact 
positioning  and  separation  of  the  energy  levels,  and  in  the  shapes  of  the  potentials  of 
the  ^Pi/2,3/2  states  which  are  involved  in  the  collisional  relaxation  process. 

This  situation  is  not  unique  to  DPALs.  In  many  systems,  one  often  arrives  at 
a  decision  point  where  several  atoms  in  one  column  of  the  periodic  table  nearly  all 
appear  equally  suited  for  a  particular  task.  It  is  often  difficult  to  forecast  which  atom 
will  be  most  effective,  but  it  is  difficult  to  switch  atoms  in  the  middle  of  the  exper¬ 
iment,  often  requiring  entirely  new  systems  and  significant  effort  fine-tuning  these 
systems.  Although  there  are  generally  several  factors  which  influence  the  decision, 
often  the  collisional  parameters  are  important.  The  theory  of  collisional  cross  sec¬ 
tions  has  been  an  area  of  intense  research  over  the  past  several  decades,  with  varying 
degrees  of  success.  In  systems  where  collisional  effects  play  a  substantial  role,  the 
method  in  this  thesis  may  be  used  to  obtain  collisional  cross  sections  between  two 
states,  given  the  potentials  for  each  state  and  the  interaction  potential.  It  is  hopeful 
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that  determination  of  collisional  cross  sections  given  input  potentials  of  a  system  may 
provide  a  partial  answer  to  this  question  in  systems  where  collisions  play  a  major 
role,  such  as  with  DPALs  [14],  or  with  most  chemical  physics  interactions. 

1.2  Background 

In  order  to  calculate  collisional  cross  sections  between  two  coupled  states,  the 
Channel  Packet  Method  is  used  [18]  [22],  The  method  employed  is  Gaussian  wave 
packet  propagation  through  a  diabadic  surface  using  the  Channel  Packet  Method [4] 
to  compute  S-matrix  elements.  This  involves  propagating  a  Gaussian  wave  packet 
forward  and  backward  in  time  via  the  Mpller  operators,  then  observing  the  overlap, 
which  is  proportional  to  the  scattering  matrix  element.  These  scattering  matrix 
elements  are  then  used  to  determine  the  collisional  cross  sections  for  several  energies. 

1.2.1  Propagation. 

As  discussed  in  Chapter  II,  the  Channel  Packet  Method  is  used  to  propagate  a 
wave  packet.  At  each  time  step,  there  are  2  Fourier  transforms  required  per  surface 
(4  total  for  a  2-state  problem).  At  each  time  step,  the  wave  packet  is  projected  on 
to  the  Mpller  hnal  state.  The  initial  and  hnal  wave  packets  are  judiciously  chosen  to 
provide  the  collisional  cross  sections  within  an  energy  range  of  interest. 

1.2.2  Collisions. 

One  can  represent  the  interaction  between  states  as  an  operator  S,  known  as  the 
scattering  operator.  This  operator  is  not  time-dependent  because  it  is  dehned  as 
taking  reactants  from  the  inhnite  past  (|t/’m))  and  resulting  in  products  in  the  inhnite 
future  il'ipout))- 
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I'lpf)  =  Sl'ipi) 


(1) 


From  this  equation,  it  is  obvious  that  interactions  between  potentials  must  be 
contained  in  S.  If  it  is  assumed  that  only  inelastic  collisions  between  the  two  nearest 
states  can  cause  an  interaction  between  the  potentials,  the  problem  of  calculating 
the  collisional  transition  rate  between  two  states  is  eqnivalent  to  the  problem  of 
calculating  the  matrix  elements  of  S'  in  a  certain  basis  [6]  [21], 

In  the  simplification  of  allowing  only  interactions  between  two  neighboring  states 
and  a  single  k  value  (plane  wave),  the  S-matrix  can  be  understood  as  a  2x2  matrix. 
(The  S-matrix  is  simply  the  representation  (/els'! A;) [6].)  Under  these  conditions,  the 
S-matrix  can  be  written  as: 

S+k,+k  S+k,-k 

^—k,+k  ^—k,—k 

Which  wave  function  is  chosen  as  the  initial  Mpller  state,  and  which  wave  function 
is  the  final  Mpller  state?  It  is  actnally  np  to  you  to  decide,  within  some  rather  flexible 
bonndaries.  The  choice  will  determine  the  ontpntted  energies  that  are  not  overrun 
by  numerical  error,  and  which  scattering  matrix  element  is  calculated.  This  arbitrary 
choice  of  wave  function  is  possible  because  the  influence  of  the  expansion  coefficients 
for  each  k  valne  are  removed  from  the  S-matrix  element  provided  the  initial  Mpller 
state  is  either  entirely  in  +k  space  or  in  —k  space,  as  discussed  in  Chapter  III. 
Similarly,  the  final  Mpller  state  must  be  either  entirely  in  +k  space  or  in  —k  space. 
For  a  spherically  symmetric  problem,  the  initial  Mpller  state  is  often  in  —k  space  and 
the  final  Mpller  state  is  often  in  +k  space.  If  the  initial  and  final  Mpller  states  are 
placed  in  the  asymptotic  limit  of  the  potential,  Gaussian  wave  packets  may  be  used. 
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Table  1.  Summary  of  Atomic  Units 


Unit  Type 

Atomic  Units 

SI  Units 

Length 

1  Bohr 

5.291772108(18)  x  10^^  m 

Mass 

1  au  (of  mass) 

9.1093826(16)  x  lO^^  kg 

Time 

1  au  (of  time) 

2.418884326505(16)  x  10^^  s 

h 

1  au  (of  angular  momentum) 

1.05457168(18)  x  lO^^  J  s 

Energy 

1  Hartree 

4.35974417(75)  x  10^*  J 

Cross  Section 

1  Bohr^ 

0.2800285204(49)  x  lO'^®  cm^ 

1.2.3  Units. 

Atomic  units  are  used  in  this  thesis  unless  otherwise  noted.  Conversion  from 
atomic  units  to  SI  units  are  summarized  in  Table  1. 

1.3  Problems  and  Objectives 

It  was  hrst  necessary  to  code  a  Gaussian  wave  packet  propagator.  Then  the  cor¬ 
relation  function  resutls  from  the  propagator  were  used  to  implement  the  Channel 
Packet  Method  and  obtain  scattering  matrix  elements.  Through  the  Method  of  Par¬ 
tial  Waves[7],  the  procedure  was  repeated  for  several  different  effective  potentials. 
The  resulting  S-matrix  elements  were  snmmed  to  obtain  the  collisional  cross  section 
for  each  energy  with  numerically  stable  data.  The  code  was  verihed  under  several 
conditions  where  the  solutions  to  the  scattering  problem  are  known.  This  code  was 
execnted  using  adiabatic  potentials  for  LiHe  from  another  Ph.D.  stndent  at  AFIT  [5]. 
From  these  calculations,  the  following  observations  were  made: 

•  The  most  difficult  hurdle  to  overcome  in  this  process  is  determination  of  the 
conpling  potential  between  the  two  states  of  the  system 

•  For  these  systems,  the  effective  potential  died  off  very  slowly  for  /  ^  Imax- 
This  reqnired  placing  the  wave  packet  far  from  the  potential  and  truncating  the 
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effective  potential  at  a  point  when  it  was  relatively  small  compared  to  other 
energy  scales  in  the  problem.  This  preserves  the  approximation  of  the  Mpller 
states  as  Gaussian  potentials. 

1.4  Scope  and  Limitations  of  Research 

The  work  performed  assumes  the  potentials  for  the  system  are  known.  These 
potentials  are  assumed  to  be  symmetric  central  potentials,  with  the  potentials  passed 
as  input  to  the  system.  Minimal  effort  was  expended  in  understanding  how  the 
parameters  for  the  potentials  are  obtained  for  a  specihc  system.  Furthermore,  it  is 
assumed  that  only  2  states  influence  the  collisional  transition  process. 

Initially,  line  broadening  was  also  to  be  investigated  using  the  Channel  Packet 
Method  and  a  review  paper  [3] ,  but  this  was  abandoned  due  to  lack  of  time.  A  study 
of  line  broadening  using  the  Channel  Packet  Method  could  be  benehcial  to  DPAL 
work  and  the  general  spectroscopy  community. 

1.5  Significance  of  Research 

Often,  one  reaches  a  point  where  several  atoms  or  molecules  appear  to  be  suitable 
for  a  given  purpose,  and  it  is  desirable  to  know  in  advance  which  atom  or  molecule 
would  be  best  suited  for  the  task  at  hand. 

In  DPALs,  one  of  several  factors  to  consider  when  selecting  an  alkali  is  maximizing 
the  ^^3/2  — ^  ^-Pl/2  cross  section  in  order  to  more  easily  obtain  a  population  inversion 
in  the  ^Pi/2  state  [14].  It  is  hoped  that  a  more  rigorous  understanding  of  the  methods 
for  computing  collisional  cross  sections  for  this  state  transition  could  better  model 
this  essential  collisional  dynamic. 

Other  applications  involving  collisional  cross  sections  arise  in  the  realms  of  chem¬ 
ical  physics  and  spectroscopy,  where  collisions  play  a  significant  role. 
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II.  Gaussian  Wave  Packet  Propagation 


Scattering  theory  describes  the  interaction  of  an  incoming  wave  packet  with  a 
potential  over  a  certain  time.  Therefore,  any  complete  discussion  of  scattering  must 
begin  with  a  discussion  of  propagating  a  wave  packet  through  an  arbitrary  potential. 
For  sake  of  convenience,  the  wave  packets  used  as  input  for  this  thesis  are  Gaussian, 
but  this  discussion  is  kept  generic  to  any  wave  function  and  any  potential. 

2.1  Theory  of  Wave  Packet  Propagation 

Central  to  any  time  evolution  problem  in  quantum  mechanics  is  the  time  evolution 
operator,  the  Hamiltonian  H\ 


H  =  f  +  V  (3) 

T  is  the  kinetic  energy  and  V  is  the  potential  energy.  In  matrix  mechanics,  the 
Schrodinger  equation  is  as  follows: 


H\i:)  =  (4) 

This  equation  has  the  following  solution  if  the  Hamiltonian  is  time-independent: 


|^(t))  =  e-'^'/^|^(0))  =  (5) 

In  the  coordinate  representation,  this  equation  has  the  following  form: 


{x\'ip(t))  =  =  {x\e  (6) 

Now,  a  problem  is  reached.  T  -|-  H  is  not  generally  diagonal  in  the  position  basis. 
Furthermore,  since  T  and  V  do  not  necessarily  commute,  the  naive  approximation 
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below  is  only  accurate  to  0(t): 


^-i{f+V)t/h  _  ^-iftlh^-iVtlh 


(7) 


A  much  better  approximation  to  make  will  be  discussed  later,  but  it  is  not  central 
to  the  discussion  at  present;  in  order  to  more  clearly  illustrate  other  features  of  the 
problem,  I  will  use  this  approximation  and  later  discuss  how  to  make  a  more  robust 
approximation.  This  approximation  has  only  solved  half  of  the  problem.  Although  it 
is  true  that  V  is  usually  diagonal  in  position  space,  T  is  not.  Fortunately,  there  is  an 
easy  way  around  this  problem.  The  kinetic  energy  of  a  particle  is  as  follows: 


2m  2m 


(8) 


This  means  T  is  diagonal  in  the  momentum  (p)  representation,  or  analogously, 
in  /c-space.  If  the  potential  is  taken  to  be  diagonal  in  the  position  representation 
(e.g.,  V{x)  is  specihed),  and  if  the  transform  from  the  x  representation  to  the  k 
representation  is  known,  then  the  problem  will  be  solved.  By  inserting  completeness 
in  the  form  of  |a;)(a;|,  one  can  see  the  transformation  of  ip  from  position  space  to 
/c-space  is  simply  the  Fourier  transform  [8]  [7] : 


dx  = —^=  f  e  ^^^ip{x)  dx  =  J^{ip{x))  (9) 

v2n  J 

The  transformation  of  tp  from  /c-space  to  position  space  is  derived  analogously, 
and  is  the  inverse  Fourier  transform  This  differs  from  the  Fourier  Transform  JF 
only  by  replacing  the  term  of  under  the  integral  with  the  term  Therefore, 
we  now  have  a  mechanism  to  evaluate  Equation  (6)  through  the  use  of  completeness 
over  k  and  x  twice  each: 


ip{k)  =  {k\ip)  =  /  {k\x){x\ 
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Wit))  =  {x 


\k){kW^tW{k'W){x'W'^W)W'\i’iO))dk  dk'dx'dx'‘ 


{x\k)e  *^‘'^^*'^^(A;|/c')(/c'|a;')e  \x”)  (x”  \^jJiO))  dkdk'  dx  dx‘ 


{x\k)e-^^^^^^/^{k\x')e-^^^^'^^l^{x'\W))  dk  dx' 


(10) 


2.1.1  Split  Operator  Method. 

As  pointed  out  previously,  it  is  necessary  to  make  an  approximation  to 
as  well  for  a  general  potential.  Up  to  this  point,  we  have  used  the  approximation  in 
Equation  (7)  for  simplicity,  despite  its  numerical  instability.  It  is  now  time  to  £x  that 
simplification.  A  much  better  approximation  is  as  follows: 


^-i{f+v)t/h  _  ^-ivtl{2h)^-iftlh^-ivtl(2h) 


(11) 


This  is  known  as  the  Split  Operator  method  [4],  which  is  widely  used  for  Time- 
Dependent  problems  because  of  its  increased  numerical  stability  at  relatively  low  cost. 
As  a  result.  Equation  (5)  is  rewritten  as: 


It  can  be  shown  that  the  equation  is  approximately  true  to  O(t^)  (that  is,  the  error 
term  scales  as  t^)  by  simply  taking  a  Taylor  expansion  of  both  sides  of  the  equation 
and  comparing  powers  of  f,  taking  care  to  consider  T  and  V  do  not,  in  general, 
commute.  This  introduces  an  algorithmic  error  of  0{t^).  This  modihcation  also  does 
not  fundamentally  alter  the  development  undertaken  in  the  previous  section.  It  does 


add  another  elementwise  multiply  (which  takes  time  0{N))  by  modifying  Equation 
(10)  to  become: 


ij{x,  t)  =  Q^j  I 

Since  the  EFT  routine  is  0{N  log2  N),  adding  an  0{N)  operation  does  not  signif¬ 
icantly  alter  runtime,  and  the  added  numerical  stability  ensures  rapid  convergence, 
so  the  use  of  this  approximation  instead  of  the  more  straightforward  Equation  (10) 
is  justihed. 

2.1.2  Absorbing  Boundary  Conditions. 

For  potentials  of  interest,  it  may  be  necessary  to  propagate  the  wave  packet  over 
long  time  scales  in  order  to  fully  pass  through  the  potential.  Through  standard 
propagation  techniques,  this  would  require  a  large  enough  grid  in  position  space  to 
hold  the  entire  packet  through  the  calculation.  However,  it  will  be  shown  in  Chapter 
III,  Equation  (25),  that  the  wave  packet  information  is  no  longer  required  once  it  has 
propagated  through  the  potential  and  beyond  a  certain  point. 

For  this  reason,  without  any  loss  of  information,  it  is  possible  to  prematurely 
truncate  the  grid  in  position  space  by  utilizing  absorbing  boundary  conditions.  The 
simplest  way  to  do  this  is  by  adding  an  absorbing  potential,  such  that  the  potential 
is  of  this  form: 


V{x) 


j  Vo{x)  -1-  iVabs{x)  for  a;  >  0 

I  Vo{x)  -  iVabs(x)  for  X  <0 


(14) 


A  natural  choice  of  the  absorbing  potential  is  a  Gaussian.  It  is  possible  for  this  ef¬ 
fective  potential  to  cause  reflections.  These  reflections  may  be  minimized  by  choosing 
a  relatively  broad  Gaussian  in  position  space.  The  effect  of  this  absorbing  potential  is 
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to  rapidly  reduce  the  value  of  the  wave  packet  at  the  boundary,  as  shown  by  Equation 

(13). 


2.1.3  3-D  Spherical  Central  Potentials. 

When  propagating  through  a  spherically  symmetric  central  potential  and  solving 
Schroedinger’s  Equation,  one  can  assume  the  following  form  for  the  solution  (where 
r  is  a  3-D  vector  and  r,  6,  0  are  the  three  scalar  components  of  that  vector): 

=  RnAr)Yr{0,  0)  =  -Un,i{r)Yr{9, 4>)  (15) 

r 

The  Yi^{d,  (p)  function  is  known,  and  therefore,  the  problem  has  been  reduced  from 
three  dimensions  to  one  dimension,  determining  the  value  of  Unj{r).  Any  spherically 
symmetric  central  potential  may  be  reduced  in  this  manner.  The  solution  for  Un,i{r) 
may  be  obtained  by  plugging  the  assumed  form  above  into  the  Schroedinger  equation, 
obtaining  [7]: 


/(/  +  1)^2 


V{r) 


Un,i{r)  = 


(16) 


2m  dr"^  2mr‘^ 

The  potential  V{r)  must  decay  faster  than  1/r  for  this  method  to  succeed.  This 
means  the  Coulomb  potential  may  not  be  analyzed  using  this  split  operator  method 
in  a  single  dimension. 


2.2  Initial  Wave  Packet 

There  is  a  remaining  issue  with  wave  packet  propagation  yet  to  be  discussed.  What 
is  |'^(a;,0))?  In  this  thesis,  Gaussian  wave  packets  will  be  used,  with  the  following 
form: 
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^  [  4{aJJ  ■  “^“*1 

The  parameters  of  the  Gaussian  wave  packet  are  as  follows: 

•  Xo  is  the  location  of  the  peak  in  position  space 

•  Axo  is  the  width 

•  ko  is  the  location  of  the  peak  momentum  in  /c-space 

•  m  is  the  mass 

The  width  of  the  Gaussian  in  /c-space  at  t  =  0  is  such  that  AxAk  =  1/2,  the 
minimum  uncertainty,  as  can  be  obtained  analytically  from  the  Fourier  transform  of 
the  above  equation: 

■^(/c)  =  Qxp  [—[AxQ{k  —  /co)]^  —  ikxo]  (18) 

There  is  flexibility  in  the  choice  of  these  parameters  for  the  Gaussian.  (In  fact, 
the  Gaussian  shape  is  not  special;  any  initial  wave  packet  can  be  chosen,  so  long  as  it 
meets  a  few  relatively  lenient  constraints  described  in  Ghapter  III.)  In  the  Ghapter 
IV,  we  will  see  that  judicious  choice  of  these  parameters  can  ensure  the  most  relevant 
collisional  information  is  extracted.  Although  this  choice  of  \ijj)  appears  arbitrary  at 
present,  it  is  because  the  scattering  operator  is  independent  of  \ip),  as  will  be  shown 
in  Ghapter  III. 
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III.  One  Dimensional  Collisional  Transitions 


3.1  Overview 

Now  that  an  arbitrary  wave  packet  may  be  propagated  through  an  arbitrary  po¬ 
tential,  the  next  step  is  to  use  the  propagation  information  to  extract  collisional  data. 
Let  S  be  an  operator  such  that: 


1^/)  =  (19) 

Here,  \'ipf)  is  the  hnal  (product)  wave  packet  and  \'ipi)  is  the  initial  wave  packet.  S 
is  known  to  be  time-independent.  The  operator  is  also  clearly  unitary  and  Hermitian. 
The  question  of  collisions  then  is  reduced  to  a  question  of  what  the  operator  S  consists 
of  for  a  given  system. 

The  S-matrix  is  dehned  as  (/c'|^|A;)  [6],  or  the  representation  of  the  operator  S  in 
k-space.  For  propagation  through  free  space  {V  =  0),  it  is  obvious  that  no  collision 
can  possibly  occur.  This  means  S  =  1,  and  k'  =  k  for  nonzero  entries  in  the  S-matrix. 
Propagation  through  an  arbitrary  potential  can  lead  to  one  or  more  of  the  following 
for  different  energies,  and  thus  different  values  of  k; 

•  Transmission  [k'  =  k) 

•  Reflection  {k'  =  —k) 

•  Energy  Transfer  {\k'\  ^  \k\) 

Energy  transfer  can  only  occur  in  an  inelastic  collision,  by  dehnition.  Although 
inelastic  scattering  is  at  the  heart  of  this  thesis,  it  is  important  to  hrst  describe  elastic 
scattering,  both  as  a  test  to  ensure  functional  code  and  as  a  precursor  for  the  more 
complex  case  of  scattering  through  an  arbitrary  potential. 
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3.2  Elastic  Scattering 


Consider  elastic  scattering  by  a  particle  with  energy  E.  Since  the  spherical  poten¬ 
tial  has  been  rednced  to  one  dimension,  there  are  only  two  possible  resulting  values 
of  k  in  the  asymptotic  limit,  as  follows: 


±k 


^  \/2mE 

h 


(20) 


Since  there  are  only  two  possible  values  of  k,  the  complete  S-matrix  for  a  single 
initial  energy  E  may  be  represented  by  a  2x2  matrix.  For  each  energy,  the  S-matrix 
is  represented  as  a  2x2  matrix  as  follows: 


§  = 


(21) 


The  sign  represents  whether  the  k- value  is  -|-|A;|  or  —\k\.  The  first  sign  is  the 
product  (output)  k-value,  and  the  second  sign  is  the  reactant  (input)  k-value;  for 
example,  S'+_  =  (+/c|>S|  —  k).  Clearly,  transmission  is  represented  by  the  diagonal 
elements  of  this  2x2  matrix,  and  reflection  is  represented  by  the  off-diagonal  elements. 
Since  the  collision  is  elastic,  all  other  matrix  elements  corresponding  to  different  k- 
values  (and  therefore  different  energies)  must  be  0,  by  dehnition.  Therefore,  this 
2x2  matrix  must  be  unitary.  The  only  remaining  step  is  to  hnd  the  actual  S-matrix 
entries,  which  will  be  discussed  in  section  3.5. 

For  a  spherically  symmetric  central  potential,  this  matrix  may  be  further  sim- 
plihed.  If  the  packet  begins  in  the  asymptotic  limit  at  large  r,  then  a  packet  with 
initial  momentum  +hk  does  not  undergo  a  collisional  transition,  so  that  matrix  ele¬ 
ment  must  be  0.  Furthermore,  a  packet  with  hnal  momentum  in  the  asymptotic  limit 
of  —hk  has  not  yet  undergone  a  collision,  and  thus  is  not  the  hnal  momentum,  so 
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that  matrix  element  also  must  be  0.  This  means  only  the  matrix  element  S'+_  has  a 
nonzero  entry  for  a  spherically  symmetric  central  potential. 

3.3  Inelastic  Scattering 

In  the  case  of  inelastic  scattering,  for  a  given  initial  energy  E,  there  are  still  only 
two  possible  input  k  values,  +k  and  —k,  but  there  are  now  an  inhnite  number  of 
possible  output  k  values,  in  general,  corresponding  to  a  different  product  energy  Ef. 

Fortunately,  in  quantum  mechanical  transitions,  we  are  often  interested  only  in 
transitions  between  a  relatively  small  number  of  total  states.  In  this  thesis,  we  will 
assume  the  two-state  approximation  holds.  Therefore,  there  are  only  a  total  of  2 
possible  hnal  energies,  and  therefore  4  possible  hnal  k-values,  corresponding  to  non¬ 
zero  elements  of  the  S-matrix. 

This  leads  to  the  following  general  form  for  the  S-matrix  for  a  given  input  energy 
E: 


‘+kl,+kl 

S+kl,+k2 

S+kl-k2 

‘  —  kl,+kl 

S—kl,—kl 

S-kl,+k2 

S-kl-k2 

'+k2,+kl 

S-\-k2,—kl 

S+k2,+k2 

S+k2-k2 

'-k2,+kl 

S—k2,—kl 

S-k2,+k2 

S-k2-k2 

As  in  the  elastic  scattering  case,  this  matrix  must  be  unitary  in  the  two-state  ap¬ 
proximation.  The  upper  half  of  the  matrix  corresponds  to  the  elastic  scattering  case, 
and  the  lower  half  of  the  matrix  corresponds  to  the  inelastic  scattering  case,  so  this 
matrix  may  be  written  simply  as: 


Skl^kl  Sk2— >kl 

Skl^k2  Sk2-+k2 
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Obviously,  the  diagonal  entries  of  this  supermatrix  are  the  elastic  scattering  and 
the  off-diagonal  elements  are  inelastic  scattering. 

This  may  be  further  simplihed  for  spherically  symmetric  central  potentials,  using 
the  same  logic  as  in  the  preceding  section.  Under  a  spherically  symmetric  potential  in 
the  asymptotic  limit,  /ci  <  0  (the  initial  Mpller  packet  is  propagating  into  the  potential 
from  the  asymptotic  limit)  and  /c2  >  0  (the  hnal  Mpller  packet  has  completely  reflected 
off  the  potential). 

Therefore,  the  nonzero  elements  of  the  S-matrix  for  spherically  symmetric  central 
potentials  when  beginning  in  the  asymptotic  limit  are  simply: 

fci 

S+k2-kl 

For  a  spherically  symmetric  central  potential,  the  following  property  must  there¬ 
fore  be  true  since  the  S-matrix  is  unitary: 


\S+k2,  -kl  P  +  \S+ki  fell  |*S*e/ast^c|  1 


(22) 


{Seiasti^  =  S+ki-ki  and  Sinda^i^  ^  5+t2,-fci)-  Tliis  important  result  for  a  spheri- 
cally  symmetric  potential  in  the  2-state  approximation  may  at  first  seem  trivial,  but 
it  gives  an  important  mechanism  to  determine  the  error  in  the  method  if  both  the 
elastic  and  inelastic  scattering  results  are  calculated. 


3.4  Interpretation  of  S-matrix 

Each  element  of  the  S-matrix  is  dehned  as  follows  [6]: 


S 


i  _ 

kf,ki 


(23) 
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Here,  6i  is  a  real  number,  and  Ai  is  an  amplitude  less  than  or  equal  to  1.  If  all  of 
the  packet  transitions  from  ki  to  /c/,  the  amplitude  Ai  must  be  1,  since  the  S-matrix 
is  unitary.  As  the  magnitude  decreases,  the  probability  amplitude  of  transitioning  to 
+kf  decreases.  For  the  special  case  of  a  spherically  symmetric  potential  (/c*  <  0  and 
kf  >  0),  if  \ki\  =  \kf  \  (elastic  scattering),  the  magnitude  Ai  is  directly  proportional 
to  the  degree  of  inelasticity  of  the  transition.  If  the  transition  is  completely  elastic, 
Ai  =  1;  if  the  transition  is  completely  inelastic,  Ai  =  0.  Each  S-matrix  element  can 
be  thought  of  as  describing  the  rotation  of  an  incoming  wave  by  an  angle  5i  relative 
to  propagation  through  free  space. 

More  precisely,  assume  propagation  of  a  wave  through  no  potential  {V  =  0). 
There  will  be  no  scattering,  and  therefore  the  S-matrix  should  be  purely  diagonal, 
but  it  is  still  not  precisely  dehned  because  it  may  differ  by  an  overall  global  phase. 
Let  =  0  in  this  case  be  the  reference  for  the  angle  to  dehne  the  global  phase  factor. 
For  an  arbitrary  reflective  potential,  the  scattering  rotation  is  given  by  5i.  The  factor 
of  2  simply  takes  into  account  a  particle’s  reflection.  The  particle  traverses  the  same 
distance  in  space  twice,  once  as  an  incoming  wave  and  once  as  a  reflected  wave. 

There  are  other  definitions  of  the  global  phase,  resulting  in  other  matrices,  such 
as  the  T-matrix  and  K-matrix.  These  matrices  are  all  related  to  one  another,  and 
provide  no  additional  fundamental  information  [6].  Therefore,  the  S-matrix  is  used 
throughout  this  thesis. 

3.5  Channel  Packet  Method 

Several  methods  have  been  devised  to  calculate  S-matrix  elements.  In  this  thesis, 
a  method  devised  by  David  J.  Tannor  and  David  E.  Weeks  in  1993  [18]  called  the 
Channel  Packet  Method  is  used.  This  method  may  be  outlined  as  follows: 

•  Choose  an  initial  wave  packet  entirely  in  either  +k  or  —k  space,  and  located 
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anywhere  in  position  space 


•  Propagate  the  initial  wave  packet  under  the  Mpller  initial  operator  to  obtain 
the  Mpller  initial  state 

•  Choose  a  hnal  wave  packet  state  entirely  in  either  +k  or  —k  space,  and  located 
anywhere  in  position  space 

•  Propagate  the  hnal  wave  packet  under  the  Mpller  hnal  operator  to  obtain  the 
Mpller  hnal  state 

Propagate  the  Mpller  initial  state  forward  in  time  under  the  full  Hamiltonian: 

•  Propagate  by  At  using  the  method  described  in  Chapter  II 

•  Build  the  correlation  function  by  recording  the  projection  on  to  the  Mpller  hnal 
state  as  a  function  of  time 

Terminate  propagation  when  future  contribution  to  correlation  function  is  negligible 

•  Take  the  Fourier  transform  of  the  correlation  function  in  energy  space 

•  Compute  the  S-matrix  element  for  each  value  of  interest  using  Equation  (25) 

The  Mpller  initial  (+)  and  hnal  (-)  operators  are  dehned  as  follows  (where  Hq  is 
the  asymptotic  Hamiltonian,  and  H  is  the  full  Hamiltonian): 

=  lim  (24) 

t— >CXD 

Physically,  the  Mpller  initial  (hnal)  operator  propagates  a  wave  packet  backward 
(forward)  in  time  under  the  asymptotic  Hamiltonian,  then  forward  (backward)  in  time 
to  time  t  =  0  under  the  full  Hamiltonian.  If  the  wave  packet  is  entirely  contained  in 
the  asymptotic  Hamiltonian,  there  is  no  net  ehect. 
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Table  2.  Expansion  coefficients  to  determine  S-matrix  elements 


S-\-k,—k 

S—k,-\-k 

S-k-k 

r]+{k)  =  0  for: 
r]-lk)  =  0  for: 

k<0 

k<0 

k>0 

k<0 

k<0 

k>0 

k>0 

k>0 

It  is  important  to  eliminate  any  influence  of  the  initial  wave  packet  choice  on  the 
observed  S-matrix  elements.  Therefore,  to  compute  the  S-matrix  elements  at  each 
k-value,  it  is  necessary  to  normalize  by  the  following  equation: 

sc,  ^  h|fc|(t/>-|ki+) 

2T[mri*_{kf)ri+{ki) 

In  this  equation,  ip_  is  the  Mqller  hnal  state  and  is  the  Mpller  initial  state.  r]_ 
is  the  expansion  coefficient  of  tjj-,  and  77+  is  similarly  the  expansion  coefficient  of  7/’+. 
is  the  Fourier  transform  of  the  correlation  function: 


{kf\S\h 


=  j  dt  (26) 

The  energy  E  is  the  total  asymptotic  energy,  including  both  potential  energy  and 
kinetic  energy.  It  is  calculating  that  is  the  essential  component  of  Equation 

(25)  in  the  Channel  Packet  Method;  the  rest  of  the  equation  is  just  normalization. 
In  order  to  more  fully  understand  this  quantity.  Equation  (26)  may  be  rewritten  as 
[18] [21]  [22]:  (here,  *  means  complex  conjugate) 

27r777 

=  -jj^[V-{-k)v+{+k)S-k,+k  +  V-{+k)v+{+k)S+k,+k+ 
r]*_{-k)ri+{-k)S-k-k  +  ri-{+k)ri+{-k)S+k-k]  (27) 

Therefore,  at  first  glance  it  appears  possible  to  obtain  a  superposition  of  S-matrix 
elements  using  this  method.  Actually,  the  r]±{±k)  are  expansion  coefficients  that  are 
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chosen  by  the  initial  wave  packet.  By  judicious  choice  of  the  Mpller  initial  and  hnal 
states,  it  is  possible  to  obtain  a  particular  S-matrix  element.  In  fact,  since  the  initial 
and  hnal  Mpller  states  are  spread  in  /c-space,  it  is  possible  to  obtain  a  small  number  of 
S-matrix  elements  using  this  technique.  Table  2  summarizes  the  appropriate  choices 
for  the  Mpller  initial  and  hnal  states  in  order  to  calculate  the  desired  type  of  S-matrix 
elements. 

3.6  Diabatic,  Adiabatic,  and  Coupling  Potentials 

In  many  instances,  inelastic  scattering  is  of  interest.  The  simplest  approach  to 
modeling  inelastic  collisions  in  quantum  mechanics  is  to  assume  only  two  states  are 
involved  in  the  collision.  All  other  contributions  are  assumed  negligible. 

In  the  two-state  approximation,  the  Hamiltonian  is  now  a  2x2  matrix  at  each 
point  in  either  position  or  momentum  space.  Unfortunately,  the  Hamiltonian  for  a 
general  system  is  not  diagonal  in  both  position  space  and  momentum  space  through 
only  the  Fourier  transform.  This  is  due  to  an  interaction  term  which  models  the 
oh-diagonal  probability  amplitude  of  transitioning  from  state  2  to  state  1  (or  vice 
versa).  This  amplitude  is  assumed  to  be  real,  and  is  known  as  the  coupling  potential. 
The  potentials  which  diagonalize  the  Hamiltonian  in  momentum  space  are  known 
as  the  diabatic  potentials,  while  the  potentials  which  diagonalize  the  Hamiltonian  in 
position  space  are  known  as  the  adiabatic  potentials. 

Through  straightforward  linear  algebra  techniques,  it  is  possible  to  obtain  the  adi¬ 
abatic  potentials  from  the  diabatic  potentials  and  the  coupling  potential  by  obtaining 
the  eigenvalues  and  eigenvectors  of  a  2x2  matrix  at  each  point  in  space.  That  is,  for 
a  given  point  in  space,  the  diabatic  potential  takes  the  following  form  [4]: 
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^ii(r)  V,2{r) 

Vu{r)  V22{r) 

This  matrix  is  symmetric  and  real,  and  therefore  the  eigenvectors  simply  form  a 
rotation  matrix. 

The  Channel  Packet  Method  is  nnchanged  by  the  two-state  approximation,  other 
than  tracking  on  which  of  the  two  potential  snrfaces  the  Mpller  initial  and  hnal 
states  are  located.  However,  when  propagating  the  wave  packet,  exponentiating  an 
operator  makes  sense  only  in  the  diagonal  representation  of  that  operator.  Therefore, 
it  is  necessary  to  convert  from  the  diabatic  to  the  adiabatic  representation  prior  to 
applying  the  potential  energy  operator  V,  and  it  is  necessary  to  convert  from  the 
adiabatic  representation  to  the  diabatic  representation  prior  to  applying  the  kinetic 
energy  operator  T  [4].  This  ensures  the  operators  T  and  V,  which  are  now  2x2 
matrices  at  each  point  in  space,  are  diagonal  prior  to  exponentiation. 

Since  both  the  adiabatic  and  diabatic  representations  are  required  to  perform  the 
wave  packet  propagation,  it  is  natural  to  ask  what  minimum  amount  of  information  is 
sufficient  in  order  to  obtain  both  the  adiabatic  and  diabatic  representations.  Any  one 
of  the  following  triplets  of  information  is  sufficient  in  order  to  obtain  the  full  diabatic 
and  adiabatic  potentials: 

•  Both  diabatic  potentials  and  the  coupling  potential 

•  Both  adiabatic  potentials  and  the  coupling  potential 

•  Both  adiabatic  potentials  and  the  rotation  angle  as  a  function  of  position 

For  a  more  detailed  discussion  of  how  to  convert  between  diabatic  and  adiabatic 
representations,  to  include  a  procedure  to  estimate  the  diabatic  potentials  from  the 
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adiabatic  potentials  and  the  crossing  point,  the  reader  is  referred  to  Faist,  Levine, 
and  Bernstein  [10]  [11]. 

3.7  Determination  of  Scattering  Cross-Sections 

Scattering  cross  sections  are  a  direct  result  of  collisional  interactions. 

As  discussed  in  Chapter  II,  the  3-D  spherically  symmetric  central  potential  prob¬ 
lem  may  be  reduced  to  one  dimension.  From  Equation  (16),  the  angular  dependence 
can  be  captured  by  introducing  an  effective  potential  term  [6]  [7] : 

K„{r)  =  V{r)  +  (28) 

The  elastic,  collisional  (absorption),  and  total  cross-sections  for  a  spherically  sym¬ 
metric  potential  (where  ki  <  0  and  kf  >  0  are  the  only  nonzero  entries)  are  calculated 
as  follows  [7]: 


CXD 

=  ^5](2(+l)|l-S‘y  (29) 

oo 

<T„,  =  <T„„  =  ^^(2/  +  l)|S},y  (30) 

1=0 

‘2'k 

^  5^(2/  +  1)(1  -  S(S‘.))  (31) 

3  1=0 

In  these  equations,  it  is  assumed  that  i,  /,  and  j  are  state  labels,  and  that  f  ^  i- 
For  the  total  cross  section,  the  summation  is  over  both  j  =  i  (elastic)  and  j  ^  i 
(inelastic)  cross  sections. 

Calculation  of  <Ji^f  (/  7^  i)  is  the  primary  goal  of  this  thesis.  Theoretically,  atot 
could  be  calculated  independently  to  determine  an  error  estimate  for  Ui^f,  but  in 
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practice,  it  was  found  that  aei  (and  therefore  atot)  required  a  summation  over  a  larger 
range  of  I  values,  and  therefore,  this  approach  was  not  used  as  an  error  estimate. 
Further  details  are  discussed  in  Chapters  IV  and  V. 

3.8  Transition  Rates  and  Cross  Sections 

The  scattering  cross  sections  in  Equations  (29-31)  assume  only  a  single  energy  is 
present  in  the  system.  In  experimental  measurements,  there  is  some  distribution  of 
energies  (generally  assumed  to  be  a  Boltzmann  distribution).  Therefore,  the  results 
from  the  Method  of  Partial  Waves  must  be  convolved  with  a  Boltzmann  disbribution 
centered  at  the  desired  temperature  of  interest  to  account  for  the  statistical  mechanics 
of  the  system. 

Fortunately,  using  the  Channel  Packet  Method,  the  S-matrix  element  for  a  wide 
range  of  energies  may  be  calculated  simultaneously.  The  data  presented  in  this  thesis 
is  the  raw  distribution  from  the  Method  of  Partial  Waves,  prior  to  convolution  with 
a  distribution  function. 

This  thesis  is  primarily  concerned  with  calculation  of  cross  sections,  from  which 
the  transition  rate  may  be  obtained  [9]. 
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IV.  Numeric  Considerations 


4.1  Overview 

The  sequential  procedure  employed  in  this  thesis  for  calculating  a  cross  section  is 
as  follows: 

•  Determine  the  range  of  energies  to  obtain  cross  section  information  across 

•  Determine  the  Hamiltonian,  wave  packet,  and  base  propagation  parameters 

•  Determine  r^ax  and  kmax  from  these  inputs  and  a  tolerable  minimum  level  of 
error 

•  Determine  Imax  using  Equation  (36) 

For  each  I  value  from  0  to  Imax- 

•  Propagate  the  wave  packet  through  this  effective  potential  and  obtain  the  scat¬ 
tering  matrix  elements  using  the  procedure  described  in  Chapter  II 

For  each  k-value  that  contains  numerically  stable  information: 

•  Combine  all  scattering  matrix  elements  to  obtain  Ui^f  using  Equation  (30) 

Continue  the  above  process  until  o'i-,f  has  converged  over  the  range  of  energies 
of  interest. 

The  remainder  of  this  chapter  focuses  on  numerical  pitfalls  and  suggestions  to 
optimize  parameter  choice,  as  well  as  a  few  suggestions  on  how  to  perform  error 
analysis. 
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4.2  General  Guidelines  on  Input  Parameters 


To  obtain  an  experimentally  observed  cross  section,  it  is  necessary  to  obtain  sev¬ 
eral  different  absorption  cross  sections  and  convolve  the  results  with  the  Boltzmann 
distribution.  Using  the  Channel  Packet  Method,  it  is  possible  even  in  the  presence  of 
numerical  error  to  calculate  many  cross  sections  from  a  single  set  of  propagation  val¬ 
ues,  provided  all  of  the  k-values  are  sufficiently  close  to  the  initial  maximum  k-value 
to  minimize  numerical  error  in  the  calculation  of  the  S-matrix  elements. 

Recall  that  AxAk  =  1/2  (atomic  units)  for  our  chosen  Gaussian  at  t  =  0.  It  is 
therefore  desirable  for  Ak  to  be  as  large  as  possible,  and  thus  Ax  as  small  as  possible, 
within  the  following  limits: 

•  The  initial  k  values  of  the  wave  packet  and  the  k  of  the  Mpller  packet  must 
either  be  all  positive  or  all  negative  (to  within  an  acceptable  degree  of  error), 
although  the  wave  packet  may  be  all  —k  and  the  Mpller  hnal  packet  all  —k,  or 
vice  versa. 

•  Smaller  Ax  means  a  larger  number  of  sampling  points  is  required  to  minimize 
propagation  error  and  obtain  an  accurate  value  for  the  projection  onto  the 
Mpller  hnal  packet. 

A  good  rough  rule  of  thumb  is  the  resulting  scattering  data  will  be  valid  in  the 
range  k^  ±  3AA;.  This  is  merely  a  rough  gideline  for  setting  parameter  information 
at  the  start,  and  is  not  guaranteed.  It  is  based  only  on  my  experience  in  analyzing 
scattering  data  in  the  absence  of  signihcant  propagation  errors,  although  there  are 
several  factors  which  may  cause  the  actual  performance  to  be  worse  (or  better). 

The  sections  that  follow  contain  analysis  to  ensure  the  code  is  functional. 
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Figure  1.  Gaussian  wave  packet  propagated  through  100  time  steps  and  compared  to 
the  analytic  solution  for  a  free  particle. 


4.3  Propagation 

To  verify  the  propagator  code  developed  operates  correctly,  two  test  potentials 
were  used  where  the  analytic  solutions  are  known  in  each  case.  Propagation  through 
free  space  {V  =  0)  is  the  simplest  case  of  wave  packet  propagation,  and  propagation 
in  the  ground  state  of  a  harmonic  oscillator  {V  =  l/2mu‘^x‘^)  results  in  a  coherent 
state. 
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Figure  2.  Wave  packet  for  diatomic  Hydrogen  propagated  to  a  time  of  1728  au  in  a 
harmonic  oscillator  as  a  coherent  state  and  compared  to  analytic  solution 


4.3.1  Free  Space  Propagation. 

In  free  space  {V  =  0),  the  time  evolution  of  the  initial  Gaussian  wave  packet  given 
by  equation  17  may  be  determined  analytically  as  follows: 


exp 


(  1  y/y  iht 

V27r(Aa:o)V  V  ^ 

-feS  +  -  xo)  - 

I  _|_  iht 
'  2(AxQ)^m 


(32) 


A  wave  packet  was  propagated  through  100  iterations  and  compared  to  this  ana¬ 
lytic  result.  Results  are  shown  in  Figure  (1). 
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4.3.2  Coherent  State  Propagation. 


To  verify  the  propagation  code  worked  with  a  potential  applied,  it  was  tested 
against  a  coherent  state.  A  coherent  state  is  a  state  in  which  the  exectation  value  is 
identical  to  the  classical  expectation  value  with  time,  both  for  the  position  and  the 
momentum,  and  the  magnitude  remains  constant  with  time  -  the  wave  packet  does  not 
spread  [8].  This  makes  a  very  convenient  test  case,  since  several  quantum  mechanical 
properties  can  now  be  compared  to  their  classical  analogue,  and  an  analytic  form  may 
be  specihed. 

To  perform  this  test,  the  following  parameters  were  used: 


Initial  Gaussian  Wave  Packet  Parameters 

xO  =  5.0000000E+000  kO  =  0 . OOOOOOOE+000 

m  =  9.1800000E+002 

<x>  =  5.0000000E+000 

dxO  =  2.3338001E-001 

Hamiltonian  Parameters 

Kinetic:  Normal 

Potential:  Simple  Harmonic  Oscillator 
m  =  9.1800000E+002  w  =  1 . OOOOOOOE-002 

Angular  Momentum:  1  =  0 

Boundary  Conditions:  None 

xO  =  O.OOOOOOOE+000 

Propagation  Parameters 

X  =-1.0000000E+002  to  1 . OOOOOOOE+002 ;  dx  = 
k  =-1.2863252E+002  to  1 . 2866393E+002 ;  dk  = 
t  =  O.OOOOOOOE+000  to  2 . OOOOOOOE+003 ;  dt  = 
(iterations:  2000) 

2.4417043E-002  (N  =  8192) 
3.1412092E-002 

l.OOOOOOOE+000 

N  is  the  number  of  points  in  the  discretization.  All  values  are  in  atomic  units. 
This  wave  packet  was  propagated  through  time  and  compared  to  the  analytic 
solution  at  a  few  data  points.  The  results  are  plotted  in  Figure  2,  where  it  is  clear 
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Figure  3.  Expectation  value  of  position  of  wave  packet  every  10  au  for  diatomic  Hy¬ 
drogen  in  a  harmonic  oscillator  as  a  coherent  state  and  compared  to  analytic  solution 


Figure  4.  Maximum  magnitude  of  wave  packet  for  diatomic  Hydrogen  in  a  harmonic 
oscillator  as  a  coherent  state.  The  value  remains  relatively  constant  around  1.307. 


the  analytic  and  computed  points  match  up  precisely.  The  expectation  value  of 
the  position  is  also  examined  as  a  function  of  time  and  plotted  in  Figure  3.  The 
expectation  value  hts  to  y  =  5cos(0.01t),  as  expected  from  the  inputs,  with  an  RMS 
error  of  4.206  x  10“^°.  Finally,  for  a  coherent  state,  the  packet  should  have  no  spread 
in  time.  Since  (t/’lt/’)  =  1  and  the  wave  packet  maintains  a  Gaussian  shape,  it  is 
sufficient  to  examine  the  maximum  value  of  the  magnitude  of  the  wave  function  as 
a  function  of  time  to  ensure  it  remains  constant.  These  results  are  plotted  in  Figure 
4.  All  of  these  tests  performed  as  expected,  indicating  the  propagation  code  worked 
properly. 

4.4  Scattering  Calculations 

To  verify  the  calculation  of  scattering  matrix  elements  after  propagation,  there 
are  two  approaches  that  may  be  used.  The  hrst  approach  is  to  simply  propagate  a 
wave  packet  through  no  potential  and  calculate  the  S-matrix  elements,  which  should 
all  be  1,  since  the  S-matrix  global  phase  factor  is  defined  by  this  value.  The  second 
method  is  to  propagate  a  wave  packet  through  a  square  potential  without  absorbing 
boundary  conditions  until  the  packet  has  propagated  completely  through  the  potential 
and  ensuring  the  square  of  the  S-matrix  elements  match  with  the  ratio  of  transmitted 
packet  to  initial  packet,  and  reflected  packet  to  initial  packet.  Since  the  propagator 
code  was  verihed  in  the  previous  section  using  a  separate  technique,  this  shows  the 
scattering  calculations  are  performed  correctly. 

4.4.1  “Scattering”  through  free  space. 

The  following  parameters  were  used  to  calculate  the  scattering  matrix  elements 
through  free  space: 
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Initial  Gaussian  Wave  Packet  Parameters 


xO  =  9.0000000E+001  kO  =  8 . OOOOOOOE+000  dxO  =  5 . OOOOOOOE-001 
m  =  9 . 1800000E+002  (atomic  units) 

<x>  =  9.0000000E+001 


Holier  Final  State  Parameters 


xO  =  l.lOOOOOOE+002  kO  =  8 . OOOOOOOE+000  dxO  =  5 . OOOOOOOE-001 
m  =  9 . 1800000E+002  (atomic  units) 


Hamiltonian  Parameters 


Kinetic:  Normal 
Potential:  None 
Angular  Momentum:  1  =  0 

Boundary  Conditions:  None 


Propagation  Parameters 


X  =  O.OOOOOOOE+000  to  3 . OOOOOOOE+002 ;  dx  =  3 . 6625565E-002 
(N  =  8192) 

k  =-8.5755010E+001  to  8 . 5775951E+001 ;  dk  =  2 . 0941394E-002 
t  =  O.OOOOOOOE+000  to  8 . OOOOOOOE+003 ;  dt  =  1 . OOOOOOOE+000 
(iterations:  8000) 


Clearly  in  this  case  the  “scattering”  is  elastic,  so  ki  =  kj.  The  expected  result  is 
for  the  real  part  of  the  scattering  matrix  element  to  be  1  and  the  imaginary  part  of 
the  scattering  matrix  element  to  be  0. 

Figure  5  contains  the  actual  results.  Since  Ax  =  0.5  for  this  run,  Ak  =  1  because 
AxAk  =  1/2  for  a  Gaussian  at  t  =  0.  The  results  become  less  accurate  with  increasing 
distance  from  ko-  Beyond  around  4Ak,  the  data  becomes  completely  unreliable  as 
numerical  error  due  to  the  small  values  of  the  wave  packet  and  hnite  storage  of 
double-precision  numbers  overshadow  any  computation  of  the  actual  matrix  element. 
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2  4  6  8  10  12  14 


k  (au) 

Figure  5.  “Scattering”  through  no  potential  for  a  wave  packet  centered  at  kg  =  8  with 
Afc  =  1  (Aa;  =  0.5).  Only  data  within  about  4Ak  is  close  to  the  correct  value,  with 
decreasing  accuracy  with  increasing  distance  from  fcg- 


4.4.2  Scattering  through  a  Square  Potential. 

A  simple  conservation  of  energy  calcnlation  was  performed  throngh  a  sqnare  poten¬ 
tial  to  farther  verify  the  scattering  code.  When  scattering  throngh  a  sqnare  potential 
that  is  0  on  both  sides  and  has  a  depth  —  Vq,  there  are  two  possible  resnlts:  reflection 
and  transmission. 

Assnme  the  initial  wave  packet  |V’o)  is  entirely  in  +k  space,  with  initial  momentnm 
+ko,  and  completely  ontside  the  sqnare  potential.  Propagate  the  wave  packet  throngh 
the  sqnare  potential  nntil  it  has  passed  completely  throngh  the  sqnare  potential  and  is 
either  entirely  reflected  or  entirely  transmitted;  let  this  hnal  state  be  |'0/)-  Calcnlate 
the  S-matrix  elements  for  transmission  (in  +k  space  beyond  the  sqnare  well)  and 
reflection  (in  —k  space  before  the  sqnare  well);  let  St  and  represent  these  S- 
matrices.  \ST{k)\‘^  +  |S'i^(fc)p  =  1  for  all  valnes  of  k  (there  is  no  other  possible  state 
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Position  (au) 


k  (au) 


Figure  6.  Initial  conditions  for  propagation  of  wave  packet  through  square  potential 


for  the  packet). 

Another  test  may  be  performed  if  absorbing  boundary  conditions  are  not  used. 
The  following  equations  hold  for  a  square  potential  [19]: 


^  K+fc|^/)P 
li+kli’oW 

l(+fc|^o)P 


(33) 

(34) 


Due  to  numerical  error,  it  is  expected  that  these  equations  will  hold  only  for  initial 
k  sufficiently  close  to  k^.  To  test  the  scattering  code  for  transmission,  the  following 
parameters  were  used  (as  illustrated  in  Figure  6): 


Initial  Gaussian  Wave  Packet  Parameters 


xO  =  2.8500000E+002  kO  =  4 . OOOOOOOE+000  dxO  =  1 . OOOOOOOE+000 
m  =  9 . 1800000E+002  (atomic  units) 
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<x>  =  2.8500000E+002 


Holier  Final  State  Parameters 

xO  =  3.1500000E+002  kO  =  4 . OOOOOOOE+000 

m  =  9 . 1800000E+002  (atomic  units) 

dxO  =  l.OOOOOOOE+000 

Hamiltonian  Parameters 

Kinetic:  Normal 

Potential:  Square  Well 

xmin  =  2 . 9500000E+002  xmax  =  3 . 0500000E+002 

yl  =  O.OOOOOOOE+000  yu  =  0 . OOOOOOOE+000 

Angular  Momentum:  1  =  0 

Boundary  Conditions:  None 

h  =-5.0000000E-002 

Propagation  Parameters 


X  =  O.OOOOOOOE+000  to  6 . OOOOOOOE+002 ;  dx  =  3 . 6623329E-002 
(N  =  16384) 

k  =-8.5770716E+001  to  8 . 5781187E+001 ;  dk  =  1 . 0471336E-002 
t  =  O.OOOOOOOE+000  to  1 . 5000000E+004;  dt  =  1 . OOOOOOOE+000 
(iterations:  15000) 


A  second  scattering  calculation  was  performed  with  identical  parameters,  except 
with  the  Mpller  dual  state  centered  at  xq  =  280  and  ko  =  —4  to  detect  the  reflection 
from  the  square  potential.  In  the  propagation,  it  is  instructive  to  note  a  few  features 
which  illustrate  the  importance  of  scattering  off  a  potential  that  may  not  be  intuitive 
due  to  the  wavelike  nature  of  quantum  mechanics. 

In  Figure  7,  the  expectation  value  has  reached  the  edge  of  the  potential  well.  The 
shape  of  the  Gaussian  has  changed  as  expected  because  some  of  the  packet  has  reached 
this  lower  potential,  and  thus  has  gained  energy.  However,  due  to  the  wavelike  nature 
of  quantum  mechanics,  there  has  been  some  reflection  off  the  potential  change  as  well, 
resulting  in  interference  in  the  packet  even  before  the  potential  has  been  reached,  and 
nonzero  negative  k. 
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Figure  9.  Propagation  through  square  potential  at  time  of  15000  au  with  transmission 
Mpller  state 
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Figure  10.  Propagation  through  square  potential  at  time  of  15000  au  with  reflection 
Mpller  state 
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In  Figure  8,  it  is  quite  obvious  that  only  some  of  the  wave  is  being  transmitted 
and  some  of  the  wave  is  being  reflected.  Interference  is  causing  only  certain  k  values 
to  be  present.  Some  of  the  packet  is  bouncing  around  inside  the  potential,  while  some 
of  the  packet  has  been  able  to  escape  from  the  potential.  The  resulting  transmitted 
and  reflected  packets  are  taking  shape,  and  the  transmitted  packet  is  beginning  to 
reach  the  Mpller  hnal  state. 

Figures  9  and  10  show  the  propagation  at  hnal  time  with  the  transmission  and 
rehection  Mpller  hnal  states,  respectively.  Under  this  set  of  parameters,  the  trans¬ 
mission  Mpller  hnal  state  exactly  matches  the  initial  wave  packet  in  fc-space,  and  the 
rehection  Mpller  hnal  state  is  the  same  as  the  initial  except  at  —k.  The  packet  is, 
in  each  case,  beyond  the  Mpller  hnal  state  in  position  space,  and  at  no  point  does 
it  possess  more  magnitude  in  |/c|  than  it  originally  had,  illustrating  the  collision  was 
elastic.  Certain  values  of  |/c|  were  rehected,  while  other  values  of  |A;|  were  transmitted, 
depending  on  the  structure  of  the  square  potential.  Equations  (33-34)  may  now  be 
understood  visually  -  it  is  the  ratio  of  the  current  value  with  the  initial  value  that  pro¬ 
vides  the  transmission  coefficient  (and  similarly  in  —k  space,  since  it  E  =  tx^k"^ /{2m), 
so  the  sign  of  k  is  irrelevent  in  determining  the  energy). 

Figures  11  and  12  show  the  correlation  function  for  transmission  and  rehection 
inputs,  respectively.  At  each  time  step  of  the  propagation,  the  result  in  position  space 
is  projected  onto  the  Mpller  hnal  state  in  position  space,  and  this  value  is  recorded  as 
the  value  of  As  indicated  by  Equation  (26),  it  is  the  Fourier  transform 

of  this  function  in  energy  space  that  provides  the  S-matrix  elements. 

Another  interesting  problem  is  encountered  here  from  a  numeric  perspective. 
Equally  spaced  points  in  position  and  k  are  natural  choices  for  the  propagation  of 
the  wave  packet.  At  each  time  step,  two  Fourier  transforms  are  necessary  in  the 
Split  Operator  propagation  scheme  used,  as  shown  in  Equation  (13).  Therefore,  it  is 
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Figure  11.  Correlation  function  of  transmission  through  square  potential 


Figure  12.  Correlation  function  of  reflection  off  square  potential 
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ideal  to  perform  these  Fourier  transforms  as  rapidly  as  possible.  The  fastest  known 
algorithm  today  is  the  Fast  Fourier  Transform,  which  relies  on  equally-spaced  inputs. 

However,  when  taking  the  inverse  Fourier  transform  to  calculate  the  S-matrix 
elements,  this  Fourier  transform  is  in  energy  space.  Since  E  =  h‘^k‘^ /{2m),  using 
eqnally  spaced  points  in  k  means  the  Fast  Fonrier  Transform  cannot  be  applied  in 
this  calculation  through  standard  techniques.  There  are  three  possible  solutions  to 
this  problem: 

•  Eqnation  (25)  relies  only  on  the  expansion  coefficients  of  the  initial  and  hnal 
states,  which  we  have  assnmed  are  Gaussian.  The  Fast  Fourier  Transform  may 
be  applied,  and  the  necessary  expansion  coefficients  may  be  calculated  for  each 
energy  valne  using  Eqnation  (18).  Although  not  generally  applicable,  the  code 
already  makes  the  assumption  that  the  Mpller  states  are  Gaussian  wave  packets, 
and  therefore  are  placed  in  the  asymptotic  limit  of  the  potential. 

•  Since  the  Eourier  transform  of  the  correlation  fnnction  is  one-time  work  and 
the  Fourier  transforms  for  propagation  are  every  iteration  work,  it  is  possible  to 
maintain  generality  by  simply  using  the  fnll  Fourier  transform  integral  expres¬ 
sion  instead  of  the  Fast  Fourier  Transform  to  compute  the  S-matrix  elements, 
and  use  the  Fast  Fourier  Transform  for  the  mnch  more  freqnent  propagation 
calculations.  This  maintans  generality,  at  the  expense  of  some  compntational 
cycles.  However,  it  was  determined  that  the  computational  cost  was  not  that 
great,  so  this  approach  was  taken  in  the  code  in  Appendix  B. 

•  It  seems  plausible  to  modify  the  Fast  Fourier  Transform  algorithm  to  work  with 
data  points  spaced  as  ^/E,  allowing  both  generality  and  speed  to  be  preserved, 
but  the  details  of  this  approach  were  not  explored,  as  the  runtime  with  the 
fnll  Fourier  transform  for  calculating  scattering  was  negligible  compared  to  the 
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Figure  13.  Square  of  the  magnitude  of  S  matrix  elements  from  a  square  potential  for 
transmission  and  reflection,  and  the  sum,  with  fcg  =  4  and  Ak  =  0.5. 


runtime  for  propagation. 

Figure  13  shows  the  magnitudes  of  transmission  and  reflection  along  with  the 
sum.  Again,  the  data  sums  to  1  as  expected  near  the  central  k  value  of  =  4.  Here, 
/Ax  =  1,  so  /Ak  =  1/2.  The  data  is  observed  to  be  accurate  from  k  =  2.7  to  k  =  5.3. 
Outside  that  range,  the  data  deviates  noticeably  from  this  sum. 

It  is  also  possible  to  compare  for  a  square  potential  using  Equations  (33)  and  (34). 
These  results  are  plotted  in  Figures  14  and  15  for  transmission  and  reflection  calcu¬ 
lations,  respectively.  Similar  results  are  observed.  Outside  the  range  from  /c  =  2.7  to 
5.5  signihcant  deviations  begin  to  appear,  at  hrst  gradually,  then  rapidly  growing  to 
overtake  any  semblance  of  a  reasonable  estimate  of  the  actual  value. s 
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Figure  14.  compared  to  T  from  Equation  (33)  for  a  square  potential,  with  fcp  =  4 

and  Afc  =  0.5. 


Figure  15.  compared  to  R  from  Equation  (34)  for  a  square  potential,  with  fco  =  4 

and  Afc  =  0.5. 
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4.4.3  Scattering  Through  a  Spherical  Central  Potential. 


For  a  spherical  central  potential,  the  procedure  is  similar,  and  a  detailed  example 
with  coupling  will  be  explored  in  Chapter  V.  The  key  difference  is  that  scattering  is 
completely  reffective  in  this  case.  If  there  is  no  coupling,  all  scattering  matrix  elements 
should  have  a  magnitude  of  1,  but  the  scattering  angles  will  differ.  Therefore,  the 
calculation  of  scattering  angles  becomes  essential  for  understanding  scattering  off  a 
spherically  symmetric  central  potential. 

If  coupling  between  multiple  states  is  introduced,  the  magnitude  of  the  scattering 
matrix  elements  may  be  less  than  1,  representing  the  transition  from  the  starting 
state  to  another  state.  The  magnitude  of  the  S-matrix  element  quantifies  the  extent 
of  the  transition. 

4.5  Cross  Section  Calculations 

Cross  section  calculations  follow  naturally  from  the  theoretical  considerations  pre¬ 
viously  discussed  in  Chapter  III,  with  two  minor  modifications.  First,  the  sums  in 
Equations  (29-31)  are  infinite,  and  the  wave  packet  propagation  for  realistic  potentials 
is  generally  the  most  computationally  intensive  step,  so  minimizing  the  required  num¬ 
ber  of  iterations  in  this  sum  while  still  converging  to  a  reasonable  answer  is  desired. 
Cohen-Tannoudji  provided  a  method  of  truncating  the  sum  [7]. 

From  Equation  (28),  with  increasing  /,  the  point  is  reached  where  the  following 
approximation  holds: 


IT 


eff 


/(/  +  1)^2 

2mr‘^ 


(35) 


When  this  point  is  reached,  the  potential  V (r)  is  no  longer  interrogated  by  a  wave 
with  energy  k.  Let  r^ax  be  the  maximum  r-value  with  non-negligible  potential  and 
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kmax  is  the  maximum  k- value  (in  absolute  value): 


I'max  i^^max  +  i)^k 


max'  max 


(36) 


It  is  possible,  depending  on  the  exact  potentials  and  wave  packet  dynamics,  that 
the  sum  may  be  truncated  sooner,  when  the  contributions  to  the  cross-section  are 
below  a  certain  threshold  value.  However,  the  above  analysis  provides  a  good  starting 
point  to  estimate  the  cutoff  for  the  sum  for  an  arbitrary  potential.  By  comparing  the 
results  for  /  Imax,  h  is  possible  to  determine  whether  the  solution  has  converged  or 
whether  further  propagations  are  required. 

The  other  numerical  issue  in  calculationg  cross-sections  is  due  to  Equations  (29- 
31).  These  equations  are  theoretically  never  negative,  as  expected,  since  a  negative 
contribution  to  the  cross-section  makes  no  physical  sense,  but  numerically,  even  for 
values  near  /cq,  the  magnitude  of  the  S-matrix  element  may  be  slightly  greater  than 
1,  resulting  in  a  computationally  calculated  negative  cross  section  contribution.  In 
order  to  avoid  this  effect,  it  is  necessary  to  force  the  maximum  value  of  the  magnitude 
of  the  S-matrix  elements  to  be  1.  In  Equations  (29)  and  (30),  it  is  not  necessary  to 
force  the  minimum  value  of  the  magnitude  of  the  S-matrix  elements  to  be  0,  as  |S'(A;)p 
cannot  be  negative,  even  under  ill  numeric  conditions. 


4.6  Error  Analysis 

When  considering  error,  it  is  of  vital  importance  to  recognize  that  error  is  a 
strong  function  of  k  (and  therefore  of  \fE)  due  to  the  presence  of  numerical  error  in 
the  discretization  process  and  chosen  initial  wave  packet.  The  code  contains  many 
redundancies  to  help  with  error  analysis  as  a  function  of  k  when  outputting  the  results 
for  the  S-matrix  and  cross  section  calculations. 

Some  suggested  ways  to  estimate  algorithmic  error  are  as  follows: 
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•  Calculate  the  S-matrix  elements  for  all  cases  (reflection  and  transmission,  elastic 
and  inelastic)  that  are  expected  from  the  specihed  potential.  The  sum  of  the 
square  of  the  magnitude  for  all  S-matrix  elements  for  a  given  ki  should  be  1; 
deviations  give  a  reasonable  error  estimate  for  |S'(A;)p  as  a  function  of  k. 

•  From  the  S-matrix  elements,  the  total,  absorption,  and  elastic  cross  sections  are 

all  calculated.  Theoretically,  this  is  redundant,  as  atot  =  This  extra 

level  of  redundancy  in  output  can  also  be  checked  to  further  estimate  the  error 
in  the  cross-section  calculations.  It  is  important  to  note  that  Imax  Kiay  differ 
for  a  el  and  This  method  may  only  be  used  if  both  and  aei  have 

converged. 

•  In  the  code,  the  command-line  option  -n=#  may  be  used  to  specify  a  minimal 
error  tolerance;  the  default  is  le  — 10.  This  specihes  a  lower  bound  on  the  error, 
although  the  actual  error  present  is  likely  to  be  much  larger  than  this  value  and 
depends  on  multiple  parameters.  Additionally,  -M=#  may  be  used  to  specify  a 
user-dehned  value  for  Imax- 

•  Algorithmic  error  in  propagation  scales  as  0{t^).  It  is  possible  to  estimate  the 
magnitude  of  this  propagation  error  by  comparing  the  error  in  S-matrix  elements 
calculated  with  At  =  to,  then  repeating  the  calculation  with  At  =  to/2. 

•  Propagator  error  may  also  arise  from  the  spacing  in  position  space  and  /c-space. 
To  reduce  the  spacing  in  /c-space,  it  is  necessary  to  sample  over  a  larger  space 
in  X.  To  reduce  the  spacing  in  x-space,  it  is  necessary  to  increase  the  number 
of  data  points. 
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V.  Alkali-Noble  Gas  Scattering  Cross  Sections 


5.1  Overview 

In  Diode  Pumped  Alkali  Lasers  (DPALs),  collisional  de-excitation  between  the 
^P3/2  and  ^Pi/2  states  is  responsible  for  moving  atoms  from  the  upper  pump  state  to 
the  upper  lasing  state.  To  establish  a  population  inversion  necessary  for  lasing,  this 
collisional  transition  should  be  maximized. 

A  reasonable  approximation  of  the  DPAL  system  is  an  alkali  colliding  with  an 
inert  gas.  Although  Rb  and  Cs  are  the  most  common  alkalis  for  DPAL  systems 
due  to  their  large  hne  structure  splitting,  which  allows  pumping  into  the  ‘^Ps/2  state 
without  pumping  the  ^Pi/2  state,  it  is  signihcantly  more  difficult  to  develop  potentials 
for  Rb-Noble  Gas  and  Cs-Noble  Gas  systems.  Therefore,  in  this  section,  LiHe  will  be 
considered. 

Under  this  physical  model,  it  was  assumed  that  He  is  in  its  ground  state, 
and  Li  was  already  in  the  excited  ^P3/2  state.  Therefore,  P  =  0  will  be  set  to  the 
asymptotic  energy  of  the  ^P3/2  state  for  convenience.  Only  hne  structure  coupling 
between  the  ^P3/2  and  ^Pi/2  states  were  permitted.  Goupling  between  other  states 
was  neglected,  and  the  hyperhne  splitting  of  each  hne  structure  state  was  ignored. 
This  led  to  a  two-state  system. 

The  following  tasks  were  performed: 

•  Galculate  the  cross  sections  for  absorption,  elastic  collisions,  and  total  colli¬ 
sions  for  LiHe  by  calculating  S-matrix  elements  starting  on  the  ^^3/2  state  and 
propagating  through  the  LiHe  potential. 

•  Galculate  S-matrix  elements  for  the  inelastic  ^^3/2  to  ^Pi/2  collision  for  a  single 
ehective  potential  with  all  other  inputs  the  same.  Sum  with  the  S-matrix  for 
elastic  scattering  to  obtain  a  reasonable  error  estimate  as  a  function  of  k. 
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5.2  Parameters 


Estimated  reasonable  operating  temperatures  for  the  system  are  around  500  K. 
In  particular,  it  would  be  desirable  to  obtain  cross  sections  in  the  range  from  300 
to  1000  corresponding  to  energies  of  0.000950050  to  0.00316683  Hartree.  These 
results  may  then  be  convoluted  with  a  distribution  function  to  obtain  experimentally 
signihcant  results. 

The  atomic  nuclear  mass  of  Li  is  12863.2  an,  and  He  is  7349.67  an.  This  leads  to  a 
reduced  mass  of  4677.23  an.  Under  this  mass,  using  \k\  =  V2mE/h,  \kmin\  =  2.98114 
and  \kmax\  =  5.44280.  This  indicates  a  wave  packet  centered  at  /c  ~  —4.2  with  a 
width  of  Ak  ~  0.65,  and  therefore  Ax  =  0.77  ~  0.8,  is  reasonable.  The  desired 
S-matrix  elements  are  S+k,-k  -  the  wave  packet  starts  at  large  r  where  the  potential 
is  negligible,  the  propagates  inward  to  small  r  {—k  direction)  until  it  is  reflected  out 
(in  the  +k  direction).  It  must  be  remembered  that  the  initial  Gaussian  must  be 
negligible  in  +k  space  for  the  Channel  Packet  Method  to  succeed,  as  shown  in  Table 
2.  Therefore,  Ak  cannot  be  made  too  large  or  a  superposition  of  S-matrix  elements 
will  be  obtained,  as  shown  in  Equation  (27).  Too  small  of  a  Ak  is  also  bad,  as  shown 
in  Chapter  IV,  as  numeric  error  overtakes  any  computational  accuracy  beyond  a  range 
of  roughly  k^  ±  3AA;.  Finally,  it  is  desirable  to  get  as  many  S-matrix  elements  in  the 
energy  range  of  interest  with  a  single  run  of  the  Channel  Packet  Method  in  order  to 
minimize  propagation  time,  which  is  another  reason  to  make  Ak  reasonably  large, 
provided  the  overlap  at  /c  =  0  is  negligible. 

For  a  spherically  symmetric  problem,  it  is  possible  to  estimate  the  propagation 
time  required  based  upon  the  position  of  the  initial  packets  and  the  minimum  k  value 
to  obtain  data  for.  Since  p  =  hk,  it  is  possible  to  estimate  the  minimum  velocity 
as  Vmin  =  hk/m.  Assuming  the  minimum  position  value  is  at  a;  =  0  and  the  Mqller 
initial  and  wave  packets  are  at  Xp,  then  a  rough  estimate  of  the  time  required  can  be 
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obtained  by  assuming  constant  velocity: 


2x„m 
f  - — - 

^max  ~  j-  1 

hk 

In  actuality,  since  the  effective  potential  for  large  I  will  slow  down  the  wave  packet, 
a  somewhat  larger  value  of  tmax  is  generally  required,  but  this  gives  a  rough  estimate 
of  the  value  to  use. 

5.3  Potentials 

The  potentials  are  the  most  questionable  part  of  the  set  of  inputs  for  these  LiHe 
calculations,  and  therefore  the  most  likely  source  of  error.  Numerous  attempts  have 
been  made  to  calculate  the  potentials  for  the  LiHe  interaction  [15]  [16] .  However,  it  is 
critical  to  note  that  these  calculations  generally  include  only  the  adiabatic  potentials, 
and  not  the  diabatic  potentials  or  the  interaction  potential.  Therefore,  only  part  of 
the  input  set  required  is  given  in  most  papers.  It  is  necessary  to  include  a  method  of 
estimating  the  diabatic  potentials  from  the  adiabatic  potentials.  Fortunately,  there 
is  a  method  given  in  [10]  and  [11].  This  method  was  applied  to  LiHe  in  [9].  The 
diabatic  potentials  provided  in  these  sources  are  just  an  estimate,  and  depends  upon 
knowledge  of  the  crossing  point  of  the  diabatic  potentials.  For  reasonable  propagation 
parameters  it  is  expected  that  error  in  this  interaction  potential  (and  thus  the  diabatic 
potentials)  is  the  primary  source  of  error  in  the  cross  section  calculations. 

There  are  a  total  of  6  states  comprising  the  level.  These  states  are  the  ^^1/2, ±1/2, 
^-^3/2, ±1/2,  and  ^^3/2, ±3/2-  In  the  adiabatic  representation,  these  potentials  form  the 
H^ni/2,  and  H^n3/2  states,  respectively  [9].  Therefore,  to  properly  model  this 

system,  it  appears  a  6x6  matrix  is  required.  This  may  be  simplihed  into  two  block 
diagonal  3x3  matrices  by  assuming  spin  does  not  flip  as  a  result  of  the  collision.  It  is 
only  necessary  to  solve  one  of  these  3x3  matrices. 
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Figure  16.  Adiabatic  potentials  for  Li-He 


Figure  17.  Adiabatic  potentials  used  in  propagation  scheme,  zoomed  in  at  avoided 
crossing  point 
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This  may  be  further  simplified  by  neglecting  rotational  (Coriolis)  coupling  [9]. 
This  decouples  the  A^ni/2  and  74^113/2  states.  It  is  sufficient  to  consider  the  following 
matrix  of  potentials  (where  each  potential  is  a  function  of  r): 

^3/2, 3/2  1/3/2, 1/2 

1/1/2, 3/2  bl/2,1/2 

Here,  the  subscript  3/2  refers  to  the  ^T3/2,+i/2  state,  and  the  subscript  1/2  refers 
to  the  ^^1/2, +1/2  state.  V3/2,i/2  =  ld/2,3/2,  so  a  2x2  symmetric  matrix  is  obtained  as 
described  in  Chapter  II.  The  V3/2,i/2  represents  an  approximation  to  the  derivative 
coupling  terms.  The  derivative  coupling  terms  arise  from  the  breakdown  of  the  Born- 
Oppenheimer  approximation.  The  eigenvalues  of  the  2x2  matrix  are  the  adiabatic 
potentials.  The  eigenvectors  of  this  matrix  are  described  by  a  rotation  matrix  at  each 
point  in  r: 


C(r) 


cos(6'(r))  sin(6'(r)) 

—  sin(6*(r))  cos(6*(r)) 


The  adiabatic  potentials  chosen  for  use  in  this  work  were  obtained  from  a  Ph.D. 
student  at  the  Air  Force  Institute  of  Technology  and  are  pending  publication  [5].  The 
potentials  are  plotted  in  Figures  16  and  17.  The  avoided  crossing  occurs  at  17  au,  as 
illustrated  clearly  by  Figure  17.  From  these  potentials,  it  was  possible  to  calculate  an 
estimated  interaction  potential  using  the  procedure  in  [9],  [10],  and  [11].  It  cannot  be 
overemphasized  that  this  procedure  is  approximate.  In  particular,  Faist  indicates  the 
interaction  potential  is  not  very  accurate  for  r  <<  Rx,  where  Rx  is  the  crossing  point. 
Fortunately,  this  error  should  be  relatively  small,  as  the  packet  cannot  propagate  very 
far  beyond  the  crossing  point  for  low  (~  10“^  Hartree)  energy. 

The  following  equation  for  6{r-)  was  used[10] [11]: 
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Figure  18.  0(r)  for  estimating  the  diabatic  potential  from  the  adiabatic  potential 


0(r) 


{I  +  ^  tan  ^(7)  —  I  tan  ^(2a(r  —  R^)  +  7)  for  r  <  Rq 
0  for  r  >  i?o 


where 


(38) 


^0  ^  ^  tan  +  tan  ^ 

The  parameters  a  and  7  may  be  estimated  from  the  difference  potential,  a  = 
1.2122  and  7  =  —1.484  were  used,  with  R^  =  17  for  the  crossing  point.  The  potential 
near  the  crossing  was  well-behaved.  However,  the  potential  was  strongly  modified 
even  far  from  the  crossing  point.  6{r)  is  plotted  for  the  LiHe  system  in  Figure  18. 
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5.4  Results 


5.4.1  Propagation. 

To  begin  exploring  the  cross  section  for  the  LiHe  problem,  the  following  propaga¬ 
tion  parameters  were  used: 


Initial  Gaussian  Wave  Packet  Parameters 

xO  =  4.0000000E+001  kO  =-4.2000000E+000 

m  =  4.6772300E+003  (atomic  units) 

Packet  all  on  Potential  #1  initially 
<x>  =  4.0000000E+001 

dxO 

=  8.0000000E-001 

Moller  Final  State  Parameters 

xO  =  5.0000000E+001  kO  =  4.2000000E+000 

m  =  4.6772300E+003  (atomic  units) 

Moller  state  on  Potential  #1 

dxO 

=  8.0000000E-001 

Hamiltonian  Parameters 

Kinetic:  Normal 
#  of  potentials  =  2 


Potential  #1 

Potential:  Interpolation  (Linear  Spline):  . ./LiHe_dia_p32_f ull_17.dat 

Angular  Momentum:  1  =  0 

Boundary  Conditions:  Gaussian  (Absorbing) 

Mag  =  l.OOOOOOOE+000  dxO  =  1 . OOOOOOOE+001 
Spherically  Symmetric 

Potential  #2 

Potential:  Interpolation  (Linear  Spline):  . ./LiHe_dia_pl2_f ull_17.dat 

Angular  Momentum:  1  =  0 

Boundary  Conditions:  Gaussian  (Absorbing) 

Mag  =  l.OOOOOOOE+000  dxO  =  1 . OOOOOOOE+001 
Spherically  Symmetric 
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Figure  19.  Propagation  to  t  =  45000  au  for  I  =  0.  The  absolute  value  of  the  initial 
momentum  is  the  dashed  line.  The  packet  has  mixed  strongly  due  to  the  interaction 
potential. 


Coupling  Potential 


Potential:  Interpolation  (Linear  Spline):  . ./LiHe_dia_V12_f ull_17.dat 


Propagation  Parameters 

X  =  0 . OOOOOOOE+000  to 

k  =-3.2166964E+002  to 

t  =  0 . OOOOOOOE+000  to 
(iterations:  400000) 

3.2000000E+002; 

3.2168927E+002; 

4.0000000E+005; 

dx  =  9.7659230E-003  (N  =  32768) 
dk  =  1.9634355E-002 

dt  =  l.OOOOOOOE+000 
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Figure  20.  Correlation  function  for  inelastic  scattering  with  a  Gaussian  centered  at  40 
au. 


Figure  19  contains  intermediate  results  of  the  propagation  for  diabatic  potentials 
calculated  at  f  =  45,000  au  for  I  =  0  (no  effective  potential).  At  this  point  in 
the  propagation,  a  significant  amount  of  mixing  has  occurred  between  the  upper  and 
lower  states.  Figure  20  shows  the  correlation  functions  from  t  =  55,  000  to  160, 000  au, 
where  the  maximum  correlation  occurs  on  the  inelastic  state  in  all  three  cases.  This 
reinforces  the  signihcant  difference  between  the  amount  of  wave  packet  transferred 
using  each  of  the  three  methods  of  calculating  the  diabatic  potential. 

5.4.2  Scattering. 

From  the  propagation  parameters,  scattering  matrix  calculations  were  obtained 
using  the  Channel  Packet  Method. 

For  the  case  of  no  effective  potential  modihcation  (/  =  0),  the  results  are  plotted 
in  fc-space  in  Figure  21,  and  in  the  more  conventional  energy  space  in  Figure  22. 
It  is  advantageous  to  observe  the  behavior  in  both  representations  to  perform  error 
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1.5 


Figure  21.  Magnitude  of  S-matrix  elements  as  a  function  of  ki  for  I  =  0  (no  effective 
potential  modification).  The  center  value  of  the  initial  and  final  Mpller  states  was 
|fc|  =  4.2. 


E(Hartree)  xlO“^ 


Figure  22.  Magnitude  of  S-matrix  elements  as  a  function  of  E  for  I  =  0  (no  effective 
potential  modification).  The  center  value  of  the  initial  and  final  Mpller  states  was 
|fc|  =  4.2  au,  corresponding  to  E  =  0.001885  Hartree. 
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analysis.  In  /c-space,  the  points  are  equally  spaced,  and  it  is  clear  that  the  numeric 
error  becomes  signihcant  at  approximately  ko  ±  3  an.  In  energy  space,  keeping  in 
mind  E  =  [hk)'^ / {2m) ,  the  points  are  no  longer  equally  spaced,  so  errors  at  larger 
energies  (and  thus  large  |A;|)  are  somewhat  magnihed  due  to  the  relative  sparsity  of 
data  points  as  compared  to  with  small  energy  (and  thus  small  \k\). 

By  examining  the  sum  of  the  inelastic  and  elastic  cross  sections  when  /  7^  0,  it 
was  discovered  that  even  for  small  effective  potentials,  the  Gaussian  shape  was  not 
a  good  approximation  for  the  shape  of  the  Mpller  state  centered  at  40  an.  In  Figure 
23,  the  S-matrix  elements  for  inelastic  and  elastic  scattering  do  not  always  sum  to  1. 

In  order  to  discover  the  reason  for  the  failure  for  /  7^  0,  using  a  more  coarse 
calculation  (TV  =  4096  and  Xmax  =  100),  the  center  of  the  Gaussian  wave  packet  was 
moved  to  30,  35,  40,  and  50  an,  respectively,  for  /  =  5  and  /  =  10  to  observe  the 
effects  on  the  scattering  matrix  elements.  Also,  a  truncated  effective  potential  was 
investigated,  under  the  following  form: 

{Veffix)  for  X  <  Xc 

(39) 

0  for  a;  >  0 

Here,  Xc  is  the  cutoff  point  for  the  potential.  Figure  24  contains  the  results  for 
I  =  5,  and  Figure  25  contains  the  results  for  /  =  10.  In  both  cases,  it  is  observed 
that  the  sum  of  the  inelastic  and  elastic  scattering  matrix  elements  are  1  over  a  wider 
range  as  the  wave  packet  is  moved  out,  and  for  the  truncated  effective  potential  the 
results  are  exactly  as  expected. 

This  means  the  Gaussian  wave  packet  is  not  a  good  approximation  for  the  Mpller 
initial  and  hnal  states,  even  though  the  effective  potential  is  very  small  in  this  region 
(7.34e  —  6  for  /  =  10  at  a;  =  40).  There  are  three  different  ways  to  handle  this 
problem.  The  simplest  solution  is  to  move  the  wave  packet  out  further.  While  this 
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Figure  23.  S-matrix  elements  with  different  effective  potential  I  values,  with  a  Gaussian 
centered  at  40  au. 
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1  =  5 


Energy  (Hartree) 


Truncated  I  (x^  =  40) 


Figure  24.  S-matrix  elements  for  I  =  5  with  the  Gaussian  position  varied. 
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I  =  10 


Energy  (Hartree) 


Truncated  I  (x^  =  40) 


Figure  25.  S-matrix  elements  for  I  =  10  with  the  Gaussian  position  varied. 
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solution  appears  simple  on  the  surface,  it  was  discovered  that  for  the  LiHe  interaction, 
I  values  of  around  100  were  necessary  to  capture  all  of  the  absorption  cross  section 
data.  Since  the  effective  potential  at  a  given  point  in  space  goes  as  0(/^),  this  would 
require  moving  the  wave  packet  out  substantially  in  space,  resulting  in  a  substantially 
longer  propagation.  Thus,  this  solution  was  not  used. 

Another  possible  solution  is  to  develop  the  Mpller  state  for  each  effective  potential 
and  use  that  potential  for  the  initial  Mpller  state.  This  must  be  done  for  all  values 
of  /  that  are  required  to  perform  the  calculation.  This  reduces  the  propagation  time 
problem,  but  requires  more  computational  work  up  front  to  determine  the  initial 
Mpller  state. 

Finally,  a  third  solution  is  to  simply  truncate  the  effective  potential  as  in  Equation 
(39),  which  ensures  the  Gaussian  is  in  fact  a  Mqller  state,  at  the  expense  of  introducing 
some  error  from  the  modified  effective  potential.  This  error  may  be  reduced  by  moving 
the  wave  packet  farther  out,  where  the  effective  potential  is  smaller.  This  solution 
was  pursued  by  performing  a  second  propagation. 

5.4. 2.1  Truncated  Effective  Potential. 

For  the  truncated  effective  potential  propagation,  the  following  parameters  were 
used: 


Initial  Gaussian  Wave  Packet  Parameters 

xO  =  5.0000000E+002  kO  =-4.2000000E+000 

m  =  4.6772300E+003  (atomic  units) 

Packet  all  on  Potential  #1  initially 
<x>  =  5.0000000E+002 

dxO  =  7.5000000E-001 

Holier  Final  State  Parameters 

xO  =  5.0000000E+002  kO  =  4.2000000E+000  dxO  =  7 . 5000000E-001 
m  =  4.6772300E+003  (atomic  units) 
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Holler  state  on  Potential  #2 


Hamiltonian  Parameters 


Kinetic:  Normal 
#  of  potentials  =  2 

Potential  #1 

Potential:  Interpolation  (Linear  Spline):  . ./LiHe_dia_p32_f ull_17.dat 
Angular  Momentum:  1  =  0;  r_max  =  4 . 9000000E+002 

Boundary  Conditions:  Gaussian  (Absorbing) 

Mag  =  l.OOOOOOOE+000  dxO  =  7 . OOOOOOOE+000 
Spherically  Symmetric 

Potential  #2 

Potential:  Interpolation  (Linear  Spline):  . ./LiHe_dia_pl2_f ull_17.dat 
Angular  Momentum:  1  =  0;  r_max  =  4 . 9000000E+002 

Boundary  Conditions:  Gaussian  (Absorbing) 

Mag  =  l.OOOOOOOE+000  dxO  =  7 . OOOOOOOE+000 
Spherically  Symmetric 

Coupling  Potential 

Potential:  Interpolation  (Linear  Spline):  . ./LiHe_dia_V12_full_17.dat 


Propagation  Parameters 


X  =  O.OOOOOOOE+000  to  6 . OOOOOOOE+002 ;  dx  =  7 . 3251129E-002  (N  =  8192) 

k  =-4.2877505E+001  to  4 . 2887976E+001 ;  dk  =  1 . 0470697E-002 

t  =  O.OOOOOOOE+000  to  5 . OOOOOOOE+006 ;  dt  =  1 . OOOOOOOE+001 

(iterations:  500000) 


The  key  differences  are:  propagation  over  a  longer  time  with  a  larger  time  step, 
fewer  data  points  to  accommodate  the  increased  number  of  iterations,  truncation  for 
the  effective  potential  above  490  an,  and  the  wave  packet  and  Mpller  hnal  states 
centered  at  500  au.  The  truncation  ensured  the  initial  and  hnal  Mpller  states  were 
Gaussian.  Since  the  packet  was  moved  out  signihcantly,  the  effect  of  truncation  was 
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minimized.  At  490  au  with  I  =  100,  the  effective  potential  is  4.50e  —  6  Hartree. 


Energy  (Hartree) 

Figure  26.  S-matrix  element  magnitude  for  I  —  0,  truncated  effective  potential  above 
490  au,  and  a  Gaussian  centered  at  500  au. 

To  determine  the  numerically  stable  range  for  this  new  data  set,  the  sum  of  the 
magnitude  squared  of  the  scattering  matrix  elements  for  I  =  0  are  plotted  in  Figure 
26.  From  this  figure,  the  region  in  which  the  data  is  valid  to  3  significant  figures  is 
from  3e-4  to  5e-3  Hartree.  This  region  will  be  referred  to  as  the  numerically  stable 
region.  Outside  this  region,  the  data  is  erroneous  due  to  finite  machine  precision. 
Only  within  this  region  can  the  numeric  error  from  using  the  Channel  Packet  Method 
be  assumed  to  be  less  than  0.1%. 

To  conhrm  correct  operation  of  the  method,  a  comparison  was  made  for  the  sum 
of  the  elastic  and  inelastic  S-matrix  elements  at  different  I  values.  The  effect  of  I  on 
the  results  noted  in  the  previous  section  was  observed  to  completely  vanish.  Error 
was  now  entirely  due  to  numeric  considerations,  as  shown  in  Figure  27.  Notice  for 
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Figure  27.  S-matrix  elements  with  different  effective  potential  I  values,  truncated 
effective  potential  above  490  au,  and  a  Gaussian  centered  at  500  au. 
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I  =  100,  the  elastic  scattering  magnitude  is  very  close  to  1  for  all  energies  that  are 
numerically  stable,  so  this  suggests  /  =  100  is  a  good  cutoff  point  to  truncate  the 
summation  given  in  Equation  (30).  For  these  potentials,  the  lower  energies  are  the 
hrst  to  undergo  purely  elastic  scattering,  beginning  at  around  /  =  40,  followed  by 
the  higher  energies.  Since  S-matrix  elements  within  the  numerically  stable  region  all 
summed  to  1,  the  results  of  this  propagation  could  be  used  to  determine  absorption 
cross  sections. 

5.4.3  Cross  Sections  (Derivative  Coupling). 

Elastic  inelastic  (cij^/),  and  total  {<7tot)  cross  sections  were  calculated  using 

Equations  (29),  (30),  and  (31),  respectively.  To  verify  I  =  100  is  the  proper  cutoff 
point  for  the  cross  section  calculations,  the  cross  sections  were  calculated  by  setting 
Imax  in  increments  of  10  -  that  is,  Imax  =  0, 10,  20, ...,  100.  The  results  for  the  absorp¬ 
tion  cross  section  are  presented  in  Figure  28.  As  expected,  the  low  energies  converge 
for  smaller  l^ax  than  the  high  energies,  since  a  smaller  effective  potential  will  cause 
the  lower  energies  to  be  reflected. 

By  examining  the  elastic  and  total  cross  sections,  it  was  clear  that  they  had  not 
yet  converged  to  their  actual  value  for  any  energies.  Therefore,  only  absorption  cross 
sections  are  reported.  It  is  the  absorption  cross  section  which  is  of  primary  interest  in 
the  state  transition,  so  this  has  minimal  effect  on  the  results.  It  was  discovered  that 
it  is  necessary  to  go  to  much  higher  l^ax  to  obtain  elastic  and  total  cross  sections  for 
the  LiHe  potentials  specihed. 

5.5  Analysis  of  Derivative  Coupling  Approximation 

The  general  trend  for  the  collisional  absorption  cross  section  is  decreasing  with 
increasing  temperature.  Since  the  absorption  cross  section  is  a  convolution  of  this 
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Figure  28.  cn^f  calculated  from  cross  section  matrix  elements  for  various  cutoff  points. 
For  Imax  =  100,  the  cross  sections  have  converged. 
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graph  with  the  Boltzmann  distribution  and  the  resonant  peaks  are  relatively  narrow, 
it  would  be  difficult  to  focus  the  energy  only  in  one  of  the  peak  resonances.  However, 
the  collisional  cross  section  should  increase  with  decreasing  temperature  in  general. 
This  effect  is  not  so  great  for  temperatures  above  1000  K,  as  the  decrease  in  absorption 
cross  section  becomes  much  less  dramatic  beyond  that  temperature.  However,  the 
closer  to  room  temperature  the  system  is  operated,  the  larger  the  collisional  cross 
section  becomes. 

Although  these  qualitative  observations  agree  reasonably  with  literature,  at  this 
point  a  more  conventional  means  of  obtaining  the  coupling  between  two  states  was 
discovered,  and  this  method  was  disbanded  at  this  point.  The  problem  with  the 
derivative  coupling  approximation  method  used  was  that  the  derivative  coupling  is 
hypothesized  to  be  present  only  near  the  crossing  point,  but  the  approximation  yielded 
derivative  coupling  terms  far  from  the  crossing  point.  Therefore,  a  more  straightfor¬ 
ward  radial  coupling  method  was  used. 

5.6  Radial  Coupling  Between  States 

There  are  3  primary  forms  of  coupling  between  the  ^^3/2  and  ^Pi/2  states.  These 
forms  of  coupling  are:  derivative  coupling  (examined  by  approximation  in  the  previous 
section),  Coriolis  coupling,  and  radial  coupling  [9].  Coriolis  coupling  and  derivative 
coupling  are  both  neglected  in  this  section.  The  formulation  herein  draws  heavily 
from  Alexander  [2]. 

5.6.1  Theory. 

Assume  the  potentials  without  spin-orbit  splitting  are  calculated.  In  the  P-level, 
this  would  yield  2  non-unique  potentials  which  will  be  referred  to  as  Hffi  and  W^. 
The  effect  of  adding  spin-orbit  splitting  is  offset  each  of  these  potentials  by  3a/2. 
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All  of  these  variables  are  a  function  of  r  in  the  general  case,  although  the  spin-orbit 
splitting  a  may  be  approximated  as  a  constant  equal  to  the  value  in  the  asymptotic 
limit.  Ignoring  Coriolis  coupling  and  derivative  coupling,  and  examining  only  hi  >  0 
(since  hi  <  0  differs  only  by  a  negative  sign  in  the  off-diagonal  Vi2  elements),  the 
potential  is  as  follows: 

If"n  +  I  0  0 

V=  0  \{2Wn  +  Wj,)  +  %  (40) 

0  ^{Wn-W^)  \{Wn  +  2W^)-a 

It  is  important  to  note  that  coupling  in  this  case  reduces  to  simply  a  2x2  matrix 
where  only  the  electronic  angular  momentum  \  J  =  3/2,12  =  1/2)  and  \  J  =  1/2,12  = 
1/2)  couple  together,  as  follows: 

Cii  Ci2  _  \{2Wu  +  Wj,)  +  %  ^{Wu-W^) 

Vi2  V22  ^{Wn-W^)  \{Wu  +  2W^)-a 

Similar  results  (with  the  exact  same  eigenvalues)  are  observed  for  12  =  —1/2. 
In  the  radial  approximation,  the  |3/2,±3/2)  state  does  not  couple  at  all  with  the 
|l/2,±l/2)  states.  Coriolis  coupling  is  responsible  for  coupling  these  two  states, 
which  is  neglected  in  this  thesis.  Therefore,  the  cross  section  may  be  calculated  using 
this  2x2  matrix. 

In  order  to  obtain  a  cross  section  for  the  ^^3/2  — P1/2  transition,  it  is  necessary 
to  combine  the  cross  sections  in  some  manner.  This  depends  upon  the  population 
distribution  among  the  hyperhne  states.  If  the  hyperhne  distribution  were  conhned 
to  the  |3/2,±3/2)  state  and  radial  coupling  were  the  only  mechanism  present,  the 
cross  section  would  obviously  be  0.  On  the  other  hand,  if  the  hyperhne  distribution 
were  conhned  to  the  |3/2, 1/2)  state,  the  cross  section  for  the  transition  would  equal 
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the  cross  section  for  the  |3/2,l/2)  state  to  transition  to  either  of  the  |3/2,  ±1/2) 
states.  Thus,  the  question  of  hnding  the  hue  structure  cross  section  depends  upon  the 
population  distribution  among  the  hyperhne  states,  since  the  coupling  mechanisms 
between  each  hyperhne  state  (and  therefore  the  cross  sections  for  each  hyperhne  state) 
diher. 

In  the  general  scenario,  if  the  hyperhne  population  distribution  is  assumed  to  be 
P(3/2,r2),  the  following  equation  holds: 

3/2 

0'3/2^1/2(P)  =  P(3/2,  n)  [(73/2,0^1/2, 1/2(P)  ±  <73/2, 0^1/2, -1/2(P)]  (42) 

0=-3/2 

In  the  radial  coupling  approximation,  (73/2,1/2^3/2,1/2  =  o-3/2,-i/2^3/2,-i/2  and  all 
other  cross  sections  are  0,  since  there  are  no  coupling  (oh-diagonal)  terms  in  the 
potential.  Therefore,  assuming  equal  hyperhne  splitting  (P(3/2,r2)  =  1/4): 

<^3/2^1/2(-S)  =  2'^3/2,1/2^1/2,1/2(-£^)  (43) 

It  is  therefore  sufficient  in  the  radial  coupling  approximation  to  hud  o'3/2,i/2^i/2,i/2(-E), 
keeping  in  mind  that  the  approximated  cross  section  for  the  hne  structure  transition  is 
simply  half  of  that  value  if  the  population  is  distributed  equally  among  the  hyperhne 
states. 

5.6.2  Cross  Sections  (Radial  Coupling). 

The  hyperhne  cross  sections  o'3/2,i/2^i/2,i/2(-E)  were  calculated  for  Li  ±  "^He,  Li  ± 
^He,  Na  ±  "^He,  K  ±  ^He,  and  ®®Rb  ±  Ar  using  the  method  outlined  in  the  previous 
section,  except  with  the  potential  given  by  Equation  (41). 

The  results  for  Li  colliding  with  diherent  isotopes  of  He  are  plotted  in  Figure  29. 
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Figure  29.  CT3/2, 1/2^172, 1/2  for  Li  +  ^He  and  Li  +  ^He 
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Figure  30.  cr3/2,i/2^i/2,i/2  for  Li  +  ^He,  Na  +  ^He,  and  K  +  ^He 

There  is  not  a  significant  variation  with  increased  mass,  although  the  cross  section 
for  Li  +  '^He  is  approximately  1.5%  higher  than  the  cross  section  for  Li  +  ^He.  The 
potentials  do  not  vary  from  one  calculation  to  another.  This  suggests  that  cross 
sections  in  the  radial  coupling  approximation  may  be  improved  slightly  by  selecting 
the  heavier  of  two  isotopes  when  a  choice  is  available.  Further  investigation  is  required 
to  determine  whether  this  result  is  specific  to  Li  +  He  or  whether  this  is  a  general 
phenomenon. 

Results  for  Li  +  ^He,  Na  +  ^He,  and  K  +  ^He  are  plotted  in  Figure  30.  In  the  radial 
coupling  approximation,  there  is  no  significant  variance  in  the  results,  particularly  for 
Li  +  He  and  K  +  He. 

The  results  for  ®^Rb  +  Ar  are  plotted  in  Figure  31.  Here,  the  cross  section  is  signif¬ 


icantly  lower  than  the  cross  sections  calculated  for  other  isotopes,  by  approximately 
a  factor  of  10“®. 


Figure  31.  173/2, 1/2^172,1/2  for  85Rb  +  Ar 

5.6.3  Comparison  to  Experimental  Values. 

There  are  some  key  differences  between  theoretical  cross  section  calcnlations  and 
experimental  measnrement  of  cross  sections.  The  cross  sections  calcnlated  from  ex¬ 
periment  would  be  the  thermal  average  of  the  cross  sections  calculated  in  the  previous 
section,  and  therefore  would  not  be  expected  to  exhibit  rapid  oscillations.  Further¬ 
more,  assuming  the  hyperfine  population  of  the  J  =  3/2  level  are  all  equally  popu¬ 
lated,  the  experimentally  observed  cross  section  is  half  the  result  calculated  for  the 
hyperhne  o'3/2, 172^1/2, 1/2  transition,  as  described  in  Equations  (42)  and  (43). 

It  was  difficult  to  obtain  accurate  experimental  comparisons  against  values  in 
the  literature  for  the  cross  sections  for  LiHe.  One  of  the  difficulties  is  that  the  fine 
structure  splitting  of  Li  is  so  small,  so  it  compounds  the  experimental  difficulty  of 
such  a  determination.  It  is  possible  to  get  a  reasonable  estimate,  though.  In  [9],  the 
cross  section  for  the  hue  structure  transition  is  listed  as  37.3  x  10“^®  cm^  (or  37.3 
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A^)  at  556  K.  Cross  sections  for  the  ^^3/2  — ^  ^-Pi/2  collision  with  other  noble  gases 
are  given  as  falling  in  the  range  from  31.4  (LiNa)  to  66.0  A^  (LiXe)  at  this  same 
temperature. 

When  the  results  for  Li  +  ^He  in  Figure  29  are  convolved  with  a  Boltzmann 
distribution  at  556  K  and  divided  by  2  as  per  Equation  (43),  the  result  computed 
is  31  X  10“^®  cm^,  which  is  somewhat  lower  than  the  experimental  value  of  37.3  x 
10“^®  cm^,  but  relatively  close.  The  results  are  expected  to  be  somewhat  low,  since 
derivative  coupling  and  Coriolis  coupling  were  completely  ignored,  though,  so  this  is 
not  a  problem. 

For  Na  +  He,  the  calculated  results  are  again  relatively  close  to  experiment, 
although  slightly  lower  if  the  hyperhne  distribution  is  assumed  to  be  equal.  Ex¬ 
perimental  cross  sections  are  given  at  385  K  as  46.0  x  10“^®  cm^,  and  at  300  K  as 
37.7  X  10“^® cm^  [16].  The  results  calculated  in  this  thesis  yield  cross  sections  of 
21.5  X  10“^®  cm^  at  385  K  and  21.4  x  10“^®  cm^  at  400  K,  respectively.  Again,  it  is 
observed  that  the  cross  sections  are  in  the  right  ballpark,  but  are  a  bit  low,  since 
derivative  coupling  and  Coriolis  coupling  were  both  neglected. 

Finally,  experimental  results  for  this  collisional  transition  in  Rb  -|-  Ar  are  difficult 
to  calculate  precisely  since  the  transitions  are  relatively  small,  but  are  around  1  x 
10“^  X  10“^®  cm^  at  340-400  K,  as  summarized  in  [17].  This  is  again  in  the  same  order 
of  magnitude  as  the  results  from  these  calculations,  although  again  these  calculations 
are  observed  to  be  somewhat  less  due  to  neglecting  derivative  and  Coriolis  coupling. 

The  calculation  of  time  dependent  results  for  Rb  -|-  Ar  allow  one  to  explore  the 
reasons  why  it  couples  so  much  worse  than  the  Alkali  -|-  He  results.  The  diabatic 
potentials  for  Li  -|-  He  are  plotted  in  Figure  32,  and  the  diabatic  potentials  for  Rb  -|- 
Ar  are  plotted  in  Figure  33.  There  are  two  factors  which  jump  out  as  explanations 
of  this  result.  First,  the  coupling  (off  diagonal)  term  in  the  potential  is  signihcantly 
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Energy  (Hartree)  Energy  (Hartree) 


Figure  32.  Potentials  for  Li  +  He  (diabatic  representation). 


Figure  33.  Potentials  for  Rb  +  Ar  (diabatic  representation). 
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lower  in  Rb  +  Ar,  meaning  the  difference  between  Wu  and  Ws  surfaces  is  less,  and 
therefore,  according  to  Equation  41,  the  coupling  is  not  as  strong.  Additionally,  in 
the  diabatic  representation  of  Rb  +  Ar,  there  is  a  well  present  on  both  the  ^^3/2  and 
^Pi/2  states,  which  is  not  present  when  compared  to  the  Li  +  He  surfaces. 


Phase  of  S-matrix  element  (Rb  +  Ar,  1=0) 


Figure  34.  Real  part  of  the  phase  of  Rb  +  Ar  scattering  matrix  for  elastic  scattering 
and  no  Centrifugal  potential. 

To  investigate  the  possible  influence  of  the  presence  of  the  well,  the  phase  of  the 
elastic  S-matrix  element  for  I  =  0  (no  Centrifugal  potential  term)  is  plotted  in  Figure 
34.  There  appears  to  be  some  resonant  structure  associated  with  the  wells  which 
may  be  further  inhibiting  the  transition  between  the  fine  structure  states;  further 
investigation  is  warranted. 
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In  summary,  the  results  obtained  compare  favorably  to  experimental  results.  Al¬ 
though  two  forms  of  coupling  were  neglected,  the  results  were  all  within  the  proper 
order  of  magnitude,  and  were  found  to  be  close  to  experimental  values.  The  results 
were  also  shown  to  exhibit  proper  behavior  when  extending  to  Rb  -|-  Ar,  where  the 
cross  section  is  signihcantly  lower.  By  examining  the  potentials  and  the  time  depen¬ 
dent  behavior,  it  is  possible  to  explain  the  low  results  for  the  Rb  -|-  Ar  cross  sections. 
Future  work  would  include  Coriolis  and  derivative  coupling  terms  to  more  fully  model 
the  coupling  between  all  states  for  Alkali  -|-  Noble  Gas  interactions,  providing  a  more 
complete  picture  of  the  hyperhne  state  transitions. 

5.7  Error  Sources 

There  are  four  possible  sources  of  error  in  these  calculated  cross  sections,  with  the 
hypothesized  most  signihcant  source  of  error  listed  hrst  and  the  hypothesized  least 
signihcant  source  of  error  listed  last: 

•  The  inputted  adiabatic  potentials 

•  The  cutoff  of  the  effective  potential,  particularly  for  larger  Centrifugal  potentials 
(larger  /,  and  therefore  higher  energy) 

•  Numeric  error  arising  from  ensuring  |1  —  Sei{k)\‘^  >  0,  as  discussed  in  Chapter 
IV 

•  Error  in  the  discretization  process 

In  particular,  the  results  for  Li  -|-  He  were  calculated  using  much  higher  quality 
potentials  than  the  results  for  the  other  Alkali  -|-  Noble  Gas  combinations  explored. 

Since  equally  spaced  points  in  /c-space  were  used,  higher  energies  contain  fewer 
points,  so  the  lack  of  signihcant  shape  in  the  high-energy  regions  may  be  due  to 
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the  smaller  sampling  in  the  high-energy  region.  This  may  be  resolved  by  either 
using  more  data  points  or  obtaining  the  Fourier  transform  of  the  correlation  function 
(already  obtained)  at  more  energy  values  and  using  the  analytic  expression  for  a 
Gaussian  for  the  Mpller  states.  However,  the  rapid  oscillations  in  the  cross  section 
calculations  would  be  averaged  out  when  convolving  with  a  Boltzmann  distribution, 
so  the  sampling  error  is  relatively  small. 
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VI.  Conclusion 


Code  was  written  to  calculate  the  absorption  cross  sections  as  a  function  of  en¬ 
ergy  for  two  coupled  one  dimensional,  spherically  symmetric  potentials.  The  Split 
Operator  Method[4]  was  used  to  propagate  the  an  initial  Mpller  state,  assumed  to  be 
a  Gaussian  in  the  asymptotic  limit,  through  a  potential.  The  correlation  between  the 
wave  packet  and  Mpller  hnal  state,  which  was  also  assumed  to  be  a  Gaussian  when 
placed  in  the  asymptotic  limit,  was  calculated  at  each  time  step.  Using  the  method  of 
Tan  nor  and  Weeks[18],  the  correlation  function  was  used  to  obtain  scattering  matrix 
elements.  From  scattering  matrix  elements  for  several  different  effective  potential 
values  and  using  the  Method  of  Partial  Waves[7],  the  absorption  cross  section  is  cal¬ 
culated.  A  method  to  estimate  the  algorithmic  error  was  also  developed. 

6.1  General  Method 

The  procedure  developed  to  calculate  absorption  cross  sections  using  the  method 
of  Tannor  and  Weeks [18]  was  developed.  Through  judicious  choice  of  parameters,  it 
was  possible  to  obtain  cross  sections  for  several  different  energies.  The  procedure  may 
be  summed  up  as  follows: 

•  Determine  the  energies  to  obtain  scattering  cross  sections  over 

•  Determine  the  /c-values  corresponding  with  those  energies  for  the  reduced  mass 
of  the  system 

•  Position  a  Gaussian  wave  packet  spanning  those  /c- values  in  the  asymptotic  limit 
and  propagate  using  the  method  in  Ghapter  II  for  /  =  0  to  Imax 

•  Galculate  the  S-matrix  elements  for  elastic  and  inelastic  scattering 
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•  By  summing  the  magnitude  of  the  S-matrix  elements  for  elastic  and  inelastic 
scattering,  determine  the  numerically  stable  region 

•  Calculate  the  absorption  cross  section  using  Equation  (30) 

The  code  for  propagation  and  scattering  was  verihed  against  several  test  cases  to 
conhrm  correct  behavior,  as  detailed  in  Chapter  IV. 

6.2  Results 

It  was  determined  that  in  order  to  use  Gaussian  Mpller  states,  it  was  necessary  to 
either  position  the  initial  wave  packet  very  far  away  from  the  potential,  or  to  truncate 
the  effective  potential  beyond  a  certain  point.  For  LiHe,  since  the  computational  time 
required  to  propagate  for  several  different  effective  potentials  would  be  prohibitive  in 
the  time  available,  the  effective  potential  was  truncated  beyond  490  Bohr.  This 
introduced  some  error  which  could  be  eliminated  by  either  starting  the  packet  farther 
out  or  by  propagating  the  wave  packet  under  the  effective  potential  for  each  value  of 

I  from  0  to  Imax- 

Two  different  methods  of  coupling  were  examined.  Under  the  hrst  method,  cou¬ 
pling  was  estimated  from  the  spin-orbit  split  coupling.  These  results  were  analyzed 
for  the  high-quality  Li  -|-  He  surfaces,  and  produced  results  in  the  right  ballpark, 
but  the  results  were  not  extremely  satisfying  and  could  not  be  extended  to  other 
Alkali  -|-  Noble  Gas  potentials  easily,  since  the  spin-orbit  split  potentials  for  other 
combinations  were  not  readily  available. 

Using  a  method  discussed  by  Millard  Alexander [2],  it  was  possible  to  obtain  a 
more  satisfying  diabatic  potential  representation  from  the  Wu  and  Ws  surfaces.  This 
method  enabled  examination  of  Li  -|-  ^He,  Li  -|-  ^He,  Na  -|-  ^He,  K  -|-  ^He,  and 
®^Rb  -|-  Ar,  assuming  only  radial  coupling.  The  results  in  each  case  were  slightly 
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below  experimental  results  since  derivative  coupling  and  Coriolis  coupling  were  both 
neglected.  However,  in  each  case,  the  cross  section  was  close  to  experimental  results. 
Furthermore,  the  substantial  difference  between  Li  +  He  and  Rb  +  Ar  was  noticed, 
and  some  explanations  from  the  dynamical  propagation  of  the  potentials  were  noted. 

6.3  Future  Work 

There  are  several  directions  for  future  work  stemming  from  thesis  could  go,  but 
they  can  be  boiled  down  to  one  of  four  main  areas.  First,  there  is  a  need  to  use 
the  technique  applied  in  this  thesis,  using  different  atomic  systems  with  different 
potentials  and  compared  to  experimental  data.  The  primary  difficulty  is  obtaining 
reliable  potentials  for  these  systems. 

Radial  coupling  was  assumed  to  be  the  only  coupling  mechanism.  One  area  of 
future  work  would  be  to  relax  this  constraint  and  introduce  Coriolis  and/or  derivative 
coupling  terms.  This  should  yield  cross  sections  which  are  closer  to  the  experimentally 
accepted  values. 

Another  thrust  of  work  stemming  from  this  thesis  would  be  development  of  a 
library  of  Mpller  initial  and  hnal  states  which  are  propagated  out  to  inhnity  under 
the  asymptotic  Hamiltonian,  then  in  to  a  reasonable  point  (perhaps  100  Bohr)  under 
a  Hamiltonian  with  the  potential  entirely  due  to  the  effective  potential,  for  several 
different  values  of  1.  At  a  range  of  100  Bohr,  most  atomic  and  molecular  potentials 
have  reached  their  asymptotic  limits,  but  for  large  values  of  I  necessary  for  the  cross 
section  calculations,  the  effective  potential  is  still  signihcant.  Therefore,  such  a  li¬ 
brary  would  provide  a  better  starting  wave  packet  for  future  scattering  calculations 
to  be  performed,  independent  of  the  atomic  system,  and  reduce  computational  time 
required  for  one  scattering  calculation  signihcantly  -  perhaps  by  as  much  as  a  factor 
of  2  or  more.  The  savings  in  time  is  because  by  starting  closer  in  space,  the  calcula- 
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tion  does  not  have  to  extend  as  far  in  time.  The  results  would  also  be  more  accurate, 
since  a  truncation  of  the  effective  potential  would  no  longer  be  required.  This  would 
yield  more  accurate  results  than  a  truncation  technique. 

Finally,  there  are  several  inefficiencies  in  the  code  as  presently  written.  The  focus 
was  on  readable,  working  code  and  not  on  speed.  It  is  possible  that  somewhat  better 
performance  could  be  obtained  by  cleaning  up  many  of  these  inefficiencies.  Although 
the  performance  would  remain  the  same  overall  order  of  magnitude  and  is  primarily 
limited  by  the  FFT  time,  it  is  expected  that  a  performance  benefit  of  perhaps  as  much 
as  a  factor  of  1.5  to  2  could  be  obtained  by  rewriting  the  code  in  a  highly  optimized 
manner.  Some  parameters  are  calculated  that  are  not  required  for  the  scattering 
calculation  and  could  be  removed  from  optimized  code.  More  details  are  provided  in 
Appendix  B. 
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Appendix  A.  Sample  Input  Data  Files 


There  are  three  separate  data  hies  used  to  control  the  propagation,  scattering, 
and  cross  section  calculations.  The  hies  required  are  a  Gaussian  Wave  Packet  hie 
(generally  ending  in  .gauss),  a  Hamiltonian  hie  (generally  ending  in  .ham),  and  a 
Propagation  hie  (generally  ending  in  .prop).  Each  of  these  hies,  along  with  all  possible 
options  that  are  allowed,  are  discussed  below. 

In  all  data  hies,  comments  begin  with  the  !  character,  and  blank  lines  are  allowed. 
No  spaces  are  allowed  between  the  text  label  and  the  equal  sign.  All  numeric  inputs 
must  be  in  atomic  units. 

1.1  Gaussian  File 

The  Gaussian  hie  is  the  simplest  of  all  of  the  data  hies,  and  contains  4  required 
parameters.  A  sample  Gaussian  hie  is  listed  below: 

!  Position  Gaussian  Wave  Packet  on  the  left  of  the  x  grid 
x0=40 

!  Spread  in  x-space.  Note  that  dkO  =  l/(2*dx0) 
dx0=0 . 8 

!  Moving  to  the  left 
k0=-4.2 

!  Reduced  Mass  of  LiHe,  in  atomic  units 
mass=4677 . 23 

The  parameters  are  fairly  self-explanatory,  specifying  a  position,  momentum, 
mass,  and  spread  of  the  initial  Gaussian  wave  packet.  These  parameters  are  used 
with  Equation  (17)  to  specify  the  initial  wave  packet.  When  performing  a  propaga¬ 
tion  with  a  Mpller  hnal  state,  the  spread  and  mass  are  also  used  for  the  Mqller  hnal 
state  automatically.  The  Gaussian  hie  is  required  by  the  Propagator  hie. 
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1.2  Hamiltonian  File 


The  Hamiltonian  file  contains  an  option  for  specifying  whether  or  not  to  turn  on 
the  kinetic  energy,  as  well  as  several  options  for  the  potential  energy  of  either  one  or 
two  state  systems,  and  the  interaction  potential  between  these  systems.  Absorbing 
boundary  conditions  are  also  specified  in  the  Hamiltonian  hie,  as  they  are  imple¬ 
mented  as  an  imaginary  modihcation  to  the  diagonal  elements  of  the  potential.  A 
sample  Hamiltonian  hie  is  listed  below: 

!  Hamiltonian  data 

!  kinetic=none  turns  off  the  kinetic  energy,  for  testing  purposes  only. 
!  kinetic=standard  is  implied  if  there  is  no  kinetic=  line  in  the  file. 
kinetic=standard 

!  Number  of  potentials 

!  If  this  parameter  is  not  specified,  then  it  is  assumed  there  is  only 

!  one  potential. 

num_potentials=2 

!  Set  up  first  potential 

!  curr_potential  tells  the  program  which  potential  the  data  applies  to 
curr_pot  ent ial= 1 

!  potential=  specifies  the  potential  types. 

!  A  full  list  of  potential  types  is  discussed  later  in  this  section. 

!  Here,  a  linear  interpolation  between  data  points  is  specified,  with 

!  points  found  in  the  referenced  file. 

potential=interpolate 

potential.! ilename=LiHe_dia_p32_full_17 . dat 

!  Set  up  second  potential 

curr_potential=2 

potential=interpolate 

potential.! ilename=LiHe_dia_pl2_full_17 . dat 

!  Set  up  interaction  potential  between  potential  #1  and  potential  #2 
!  The  parameters  for  a  coupling  potential  are  identical  to  parameters 
!  for  an  absorbing  potential. 
coupling=interpolate 
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coupling_f ilename=LiHe_dia_V12_full_17 .dat 

!  Set  up  absorbing  boundary  conditions  (optionally) 

!  boundary=  sets  the  absorbing  boundary  condition  type  (currently, 

!  only  gaussian  or  none  is  supported) 

!  If  boundary=none  or  no  specification  is  present,  absorbing  boundary 

!  conditions  are  disabled 

boundary=gaussian 

!  It  is  possible  to  specify  a  large  absorbing  boundary  conditions  with 
!  magnitude  greater  than  1,  although  it  is  not  usually  necessary, 
boundary _magnitude=l 

!  dxO  specifies  the  width  of  the  Gaussian  absorbing  boundary  condition. 
boundary_dx0=10 


Several  different  potentials  are  supported  in  addition  to  the  linear  interpolation. 
Any  of  these  potentials  can  be  used  with  either  potential=  or  coupling=.  Examples 
of  each  potential  are  given  below.  It  is  also  possible  to  specify  potential=none 
(or  coupling=none)  for  no  potential.  When  using  two  potentials,  all  potentials  are 
specihed  in  the  diabatic  representation. 

The  Hamiltonian  hie  is  required  by  the  Propagator  hie. 

1.2.1  Simple  Harmonic  Oscillator. 

An  example  of  a  Simple  Harmonic  Oscillator  potential  specihcation  is  as  follows: 

potential=sho 
potential_x0=0 
potential_mass=0 . 5 
potential_omega=l 

These  variables  specify  a  potential  of  the  form  H  —  All  parameters 

are  required. 
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1.2.2  Square  Well. 


A  square  well  potential  may  be  specified,  with  different  heights  below,  at,  and 
above  the  well.  The  following  code  specifies  a  square  potential: 

!  Potential  energy 

potential=square_well 

!  Start  of  the  square  well 

potential_xmin=10 

!  End  of  the  square  well 

potential_xmax=20 

!  The  potential  value  inside  the  well 
potential_height=-0 . 166 

!  With  X  <  potential_xmin,  the  potential  value 
potential_lower=0 

!  With  X  >  potential_xmax,  the  potential  value 
potential_upper=0 


This  sample  sets  a  square  well  with  a  depth  of  0.166  an  inside  the  well,  and 
no  potential  outside  the  well.  All  parameters  are  required  to  specify  a  square  well 
potential. 

1.2.3  Interpolation. 

This  performs  a  linear  interpolation  between  each  point  in  the  data  file  specified. 
The  parameters  are  as  indicated  in  the  first  example  in  this  section.  There  are  two 
requirements  for  this  data  file: 

•  The  data  must  contain  two  columns.  The  first  column  is  the  x  (or  r)  value. 

The  second  column  is  the  potential  energy,  in  atomic  units. 

•  The  data  must  be  sorted  from  least  to  greatest  in  the  first  column. 

For  X  less  than  the  first  entry  in  the  file,  linear  extrapolation  is  used,  with  the 
slope  from  the  first  entry  used  to  extrapolate.  For  x  greater  than  the  last  entry  in 
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the  file,  linear  extrapolation  is  used,  with  the  assumption  of  no  slope.  (For  spherical 
potentials,  the  potential  should  level  off  at  some  value.) 

1.2.4  Morse  and  Lennard-Jones  6-12. 

Morse  and  Lennard-Jones  potentials  were  originally  planned  to  be  implemented, 
but  were  not  implemented  due  to  the  difficulty  in  htting  the  input  potentials  to  a 
Morse  or  Lennard-Jones  potential  while  maintaining  accurate  behavior.  Therefore, 
use  of  the  Morse  or  Lennard-Jones  potentials  have  not  been  fully  tested.  It  is  recom¬ 
mended  that  the  Linear  Interpolation  is  used  instead  of  the  Morse  or  Lennard-Jones 
potentials,  as  it  will  allow  a  more  arbitrary  potential  form. 


1.3  Propagation  File 


The  Propagator  hie  combines  the  initial  Gaussian  wave  packet,  under  some  Hamil¬ 
tonian,  and  discretized  under  further  parameters  as  specihed  within  the  propagation 
hie.  An  example  of  a  Propagation  hie  is  as  follows: 

!  Path  from  Propagator  file  to  initial  wave  packet 
gaussian=. . /LiHe_left .gauss 

!  Path  from  Propagator  file  to  Hamiltonian  file 
hamiltonian= . ./LiHe_2p.ham 

!  Holler  parameters  are  not  required.  If  not  specified,  there 
!  is  no  Holler  calculation,  just  a  propagation. 

!  It  is  assumed  the  spread  (dxO)  and  mass  are  the  same.  The 
!  Holler  final  state  is  assumed  to  be  a  Gaussian.  The  Holler 
!  initial  state  is  assumed  to  be  the  initial  wave  packet. 
moller_f inal_x0=50 
moller_f inal_k0=4 . 2 

!  Set  the  packet  and  Holler  final  states,  using  potential  numbers 
!  as  in  the  Hamiltonian 

!  This  is  the  elastic  scattering  case,  so  they  are  the  same 
!  These  parameters  are  not  required,  and  default  to  1  if  not  specified. 
!  Wave  packet  on  the  first  potential 
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packet _potential=l 

!  Holler  final  state  on  the  first  potential 
moller_potential=l 

!  Set  1-value  for  effective  potential  modifications. 

!  If  this  parameter  is  not  specified,  the  default  is  0. 
l_value=0 

!  Set  cutoff  position  above  which  the  effective  potential  is  not  applied. 

!  If  this  parameter  is  not  specified,  the  default  is  0 

!  (which  means  no  cutoff) . 

l_max=490 

!  Directory  to  output  resulting  calculations  to. 

!  This  directory  does  *not*  have  to  exist  prior  to  execution 
!  WARNING:  all  data  in  this  directory  will  be  erased  when  the  program 
!  is  executed  with  this  parameter  file  as  input, 
output _directory=out 

!  Other  parameters  for  propagation: 

!  Number  of  data  points  to  discretize  the  data  set  into 

num_points=32768 

!  Minimum  and  maximum  x-values 

!  NOTE:  k-value  data  is  only  accurate  when  x_min=0.  This  is  not  a  concern 

!  since  the  actual  propagations  are  spherical  potentials,  where 

!  x_min=0  obviously. 

x_min=0 

x_max=320 

!  t_min=0  is  assumed.  That  parameter  may  not  be  changed. 

!  Set  the  maximum  time  to  propagate  to 
t_max=200000 

!  Set  the  step  size  in  time 
dt=l 
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Appendix  B.  Code  Comments 


2.1  Acknowledgements 

The  author  is  grateful  to  the  developers  of  the  Fastest  Fourier  Transform  in  the 
West  library  [12],  Additionally,  the  author  is  grateful  to  the  developers  of  SCons 
version  1.2.0  [13],  which  was  used  to  automate  the  program  build  process,  and  GCC 
4.3  [1]  (gfortran  in  particular),  which  was  used  to  compile  the  program. 

2.2  Code  Availability 

The  code  is  available  upon  request  from  the  author,  who  may  be  reached  by  e-mail 
at  samuel.butler@us.af.mil,  or  from  the  author’s  thesis  advisor.  Dr.  David  Weeks. 

The  procedure  in  Chapter  IV  translates  into  the  following  for  the  user  of  the 
available  code: 

•  Determine  manually  the  required  energies,  Hamiltonian,  wave  packet,  and  base 
propagation  parameters 

•  Create  a  .gauss,  .ham,  and  .prop  file  based  on  this  information  (see  Appendix 
A). 

•  Determine  the  tolerable  lower  bound  on  the  error 

•  Run  the  following  code  to  determine  Imax  and  generate  specihc  propagator  files: 

$  prop_copy  myfile.prop  -n=le-5 

(Alternatively,  use  -yi=lmax  to  specify  a  user-dehned  maximum  1-value.)  Several 
different  propagator  files  will  be  generated  ending  with  .0  to  dmax- 

•  Propagate  each  hie  (except  the  original  master  hie)  and  calculate  S-matrix  ele¬ 
ments: 
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$  scattering  myf ile .prop . 0  -n=le-5 
$  scattering  myf ile .prop. 1  -n=le-5 

$  scattering  myf ile .prop . <l_max>  -n=le-5 

(NOTE:  By  default,  any  propagation  results  in  an  outputted  file  each  iteration. 
This  may  result  in  very  large  space  requirements  on  your  hard  drive,  and  may 
result  in  substantially  increased  runtime.  You  may  use  -x=^  or  -k=^  to  control 
this  output.  If  -x=0  or  -k=0  is  specihed,  only  the  initial  and  hnal  packets  are 
outputted.  Using  -x=1000  and  -k=1000,  for  example,  will  outpnt  x  and  k  only 
every  1000  iterations  instead  of  the  defanlt  of  every  iteration.) 

•  Run  the  following  program  to  determine  the  cross  sections: 

$  xsection  myf ile. prop  -n=le-5 

•  Finally,  use  the  cross-section  data  to  calculate  the  collisional  transfer  rate  man¬ 
ually  by  convolution  of  the  absorption  cross  section  with  the  Boltzmann  distri¬ 
bution. 

Additionally,  a  program  named  propagator  is  available.  This  program  only  per¬ 
forms  the  wave  packet  propagation,  withont  calculating  S-matrix  elements. 

2.3  Known  Issnes 

On  some  platforms,  it  was  noted  that  the  scattering  and  xsection  programs  would 
result  in  a  segmentation  fanlt  when  ontputting  certain  parameters  to  the  screen.  If 
a  problem  is  enconntered,  do  not  rnn  the  affected  program  with  the  -s  command-line 
option.  The  resnlt  of  the  program  withont  using  -s  will  appear  silent,  althongh  the 
program  will  rnn  as  expected  and  will  ontpnt  resnlts  in  the  directory  specihed  in  the 
.prop  hie. 
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In  the  code  used  for  the  propagation,  it  is  assumed  that  the  initial  wave  packet  is 
in  the  asymptotic  Hamiltonian,  outside  of  influence  of  the  potential.  Therefore,  the 
Mpller  initial  state  is  assumed  to  be  the  initial  wave  packet.  Furthermore,  this  initial 
wave  packet  is  assumed  to  be  a  Gaussian.  The  Mpller  hnal  state  is  also  assumed  to 
be  in  the  asymptotic  Hamiltonian  region,  outside  the  influence  of  the  potential,  and 
therefore  is  also  a  Gaussian. 

Gommand  line  options  are  consistent  from  one  program  to  another,  although  not 
all  programs  will  use  all  of  the  command-line  options.  This  allows  a  user  to  specify 
one  set  of  command-line  options  and  simply  pass  those  same  options  to  each  program 
being  called. 

This  is  the  author’s  hrst  major  project  in  FORTRAN,  and  as  such,  there  may 
be  some  hidden  bugs  or  other  issues  that  I  am  unaware  of.  It  is  known  that  the 
code  is  relatively  inefficient.  Some  of  the  inefficiencies  are  due  to  calculations  that 
are  nice  to  have,  but  not  essential,  in  calculating  the  cross  sections.  For  example, 
the  code  verihes  the  wave  packet  has  a  projection  of  close  to  1  at  all  iterations  in 
time,  if  absorbing  boundary  conditions  are  not  used.  This  redundancy  is  a  good 
check  to  perform  to  identify  numeric  error,  but  is  not  fundamental  to  the  calculation. 
As  another  example,  the  code  calculates  the  expectation  value  in  x  and  outputs  this 
value  to  a  hie  at  all  time  steps.  This  information  is  also  generally  not  needed,  but 
enables  some  checks  to  ensure  functioning  code  in  certain  situations.  However,  these 
calculations  do  not  affect  the  order  of  magnitude  of  the  runtime  substantially. 

If  performance  is  especially  sluggish,  one  of  two  major  factors  are  likely  to  blame. 
If  the  number  of  data  points  in  the  discretization  is  specihed  in  the  parameter  hie  as 
N,  then  for  small  N,  the  runtime  will  be  limited  by  the  I/O  of  data  hies  to  the  hard 
drive  for  ip  and  at  each  point  in  time.  To  reduce  this  strain  (and/or  the  amount 
of  disk  space  consumed),  use  the  -k=#  and  -x=pp  command  line  options  to  specify 
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outputs  only  every  7^  iterations  instead  of  the  default  of  every  iteration.  If  -k=0  or 
-x=0  is  specified,  only  the  initial  and  hnal  wave  packets  are  outputted.  (The  default 
is  -k=0  and  -x=0.) 

At  some  point,  a  majority  of  the  computational  burden  per  iteration  will  shift  to 
calculation  of  the  FFT  at  each  time  step.  To  speed  up  an  FFT-limited  routine,  use  a 
smaller  value  of  N.  If  this  is  not  acceptable,  then  either  increasing  the  time  step  or 
decreasing  the  maximum  time  will  also  enhance  run  time. 

2.4  Disclaimer 

The  code  from  the  author  is  provided  as-is,  with  no  guarantee,  implied  or  stated, 
of  working  correctly.  Although  the  author  has  much  experience  in  programming  in  C 
and  C++,  this  is  the  author’s  hrst  major  project  in  FORTRAN,  and  thus  it  is  quite 
likely  that  there  are  amateur  coding  mistakes  that  may  cause  problems  when  porting 
the  code  to  different  platforms.  The  target  operating  system  for  this  code  is  Linux; 
some  changes  are  necessary  in  order  to  run  the  code  on  another  operating  system. 
The  code  was  mostly  written  comply  with  the  FORTRAN  95  standard,  although 
some  GCC  extensions  and  FORTRAN  2003  features  may  have  also  been  used. 

The  use  of  the  code,  in  whole  or  in  part,  in  the  development  of  commercial  soft¬ 
ware  is  strictly  prohibited.  It  is  provided  strictly  for  free  personal,  academic,  and 
government  use  only. 
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